fcFOSR-TR*  8  6-  0  935 


AN  END  OF  THE  YEAR  REPORT  TO 


AIR  FORCE  OFFICE  OF  SCIENTIFIC  RESEARCH 


"  .ADQUARTERS  MILITARY  AIRLIFT  COMMAND 


for  a  Grant 

Approved  for  public  relea9i  ; 
entitled  distribution  unlimited. 


S  Id 


DEVELOPMENT  AND  EVALUATION  OF  A  CASUALTY  EVACUATION  MODEL 
FOR  A  EUROPEAN  CONFLICT 


r  5 

»**  Km  !»> 

C-<  £4  O 


«  g  «: 

M  ^ 

L  *-‘i  ♦* 

'A  -r  i 


«  t  "H  •» 


-0  «  .*>  -•  H 

-1  ,  .->  •*  -mm  . 

‘  i  ?  b 

£  r  -  0  -  q  S 


Prepared  by: 


Jeffery  L.  Kennington 
Department  of  Operations  Research 
and  Engineering  Management 


SOUTHERN  METHODIST  UNIVERSITY 
Dallas,  Texas  75275 


,*r-r*=U-'V  \ 


86  10  16 


1 4  9  ** 


December,  1985 


UNCLASSIFIED _ 

Jtcu*n>  Classification  of  this  page 


m 


c^rti  r)?>\o5 


REPORT  DOCUMENTATION  PAGE 


lb.  restrictive  MARKINGS 


2.  SECC'RITV  CLASSIFICATION  AUTHORITY 


2B  OECLASSIFICATION/OOWNGRADING  SCHEDULE 

N/A 


A  PERFORMING  ORGANIZATION  REPORT  NUMBERlS) 


3-  DiSTRl  BUT  ION/ A  V  A I  LAB  I  LIT  Y  OF  REPORT 

.Approved  for  public  release;  distribution 
unlimited 


5.  MONITORING  ORGANIZATION  REPORT  NUMBERlS. 

AFOSR-TR.  8  6-0  0  35 


6»  NAME  OF  PERFORMING  ORGANIZATION  16b.  OFFICE  SYMBO  L  7a.  NAME  OF  MONITORING  ORGANIZATION 
_  .  ,  ,  ,  ...  ,  ^  (if  applicable  i 

Southern  Methodist  University 

AFOSR  _ 


7b.  AOORESS  (City,  State  and  ZIP  Code/ 

g  Building  410 

Bolling  AFB,  D.C.  20332-6448 


a.  NAME  OF  FUNDING/SPONSORING 
ORGANIZATION 

AFOSR 


8c  ADDRESS  (City.  State  and  ZIP  Codei 

Building  410  program 

Bolling  AFB  D.C.  20332-6448  element  no. 

61102F 

II  title  (Include  Security  ClM,.«colion<DEVELOPMENT  &  EVALUA- 

CASUALTY  EVACUATION  MODEL  FOR  EUROPIAN  CONFLICT 


12.  personal  autmorisi 

ON 


8b.  OFFICE  SYMBOL  9.  PROCUREMENT  INSTRUMENT  IDENTIFICATION  NUMBER 
(If  applicable  > 

NM  AFOSR  83-0278 


10.  SOURCE  OF  FUNDING  NOS. 


PROJECT 

TASK 

WORK  UNIT 

NO. 

NO. 

NO. 

2304 

A 

13a.  TYPE  OF  REPORT 

Annual 


16  supplementary  notation 


13b.  TIME  COVERED 


14.  DATE  OF  REPORT  (Yr.,  Mo.,  Day)  I  15.  PAGE  COUNT 


>  ^0Nnv8S  December,  1985 


17 

COSATI  COOES  f 

FIELD 

GROUP 

SUB  GR.  1 

18.  SUBJECT  TERMS  (Continue  on  reverse  if  necessary  and  identify  by  block  number) 


19.  ABSTRACT  (Continue  on  reverse  i f  necessary  and  identify  by  block  numbers  _  ,  _  . 

The  mathematical  model  presented  in  this  paper  is  beyond  the  state-of-the-art  for  existing 
mathematical  programming  software.  Based  on  this  working  paper  and  discussions  with  LtCol 
McLain,  Professor  Kennington  Proposed  a  research  plan  to  refine  the  McLain  model  and  comput; 
tionally  investigate  new  algorithms  for  solving  the  model.  McLain  and  Chmielewski  collected 
the  data  and  developed  the  model  generator  while  Kennington  and  his  students  investigated 
alternative  solution  algoithms. 


20.  OISTRIBUTION/AVAILABILITY  OF  ABSTRACT 
UNCLASSIFIEO/UNLIMITEO  □  SAME  AS  RPT.  □  OTIC  USERS  □ 


121.  ABSTRACT  SECURITY  CLASSIFICATION 


22».  NAME  OF  RESPONSIBLE  INOIVIOUAL  22b  TE  LEPHONE  NUMBER 

DR  MARC  Q.  JACOBS  (Include  Arte  Code, 

am  Manager.  Mathematical  &  Information  S  :iences  767-4939 


DD  FORM  1473,  83  APR  EDITION  OF  1  JAN  73  »S  OBSOLETE. 


22c  OFFICE  SYMBOL 

NM 


SECURITY  CLASSIFICATION  OF  THIS  PAGE 


This  work  was  initiated  by  Lt.  Col.  Dennis  R.  McLain  and  Mr.  Thomas 
E.  Kowalsky  of  DCS/Plans,  Headquarters  Military  Airlift  Command  at  Scott 


AFB.  Lt.  Col.  McLain  has  transferred  to  the  Pentagon  and  the  work  at 
Scott  is  being  continued  by  Mr.  Kowalsky  and  Capt.  Robert  Chmielewski. 

The  work  originated  with  a  November  1982  working  paper  by  McLain  entitled 
"Wartime  Evacuation  of  European  Theater  Casualties  to  the  United  States". 
The  mathematical  model  presented  in  this  paper  is  beyond  the 
state-of-the-art  for  existing  mathematical  programming  software.  Based 
on  this  working  paper  and  discussions  with  Lt.  Col.  McLain,  Professor 
Kennington  proposed  a  research  plan  to  refine  the  McLain  model  and 
computationally  investigate  new  algorithms  for  solving  the  model.  McLain 
and  Chmielewski  collected  the  data  and  developed  the  model  generator 
while  Kennington  and  his  students  investigated  alternative  solution 
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This  model  is  a  member  of  the  class  of  optimization  models  known  as 
networks  with  side  constraints.  Within  this  general  area,  Kennington 
managed  three  projects  with  three  different  students.  Dr.  Ellen  Allen's 
dissertation  was  directed  toward  multicommodity  problems  and  appears  in 
Appendix  A.  A  related  study  involving  the  Equal  Flow  Problem  was 
undertaken  by  Dr.  Bala  Shetty.  This  study  is  reported  in  Appendix  B. 

New  convergence  results  related  to  both  studies  is  presented  in  Appendix 
C.  A  study  on  the  general  problem  of  networks  with  side  constraints  was 
pursued  by  Dr.  Keyvan  Farhangian.  The  results  of  that  report  are 
presented  in  Appendix  D  which  is  scheduled  to  appear  in  The  Annals  of  the 
Society  of  Logistics  Engineers. 
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During  1984,  Dr.  Narendra  Karmarkar,  of  Bell  Laboratories,  reported 
that  he  had  developed  a  new  technique  to  solve  linear  programs  based  upon 
a  projective  algorithm.  He  claimed  a  worst-case  bound  on  solution  steps 
that  was  much  better  than  the  simplex  method,  and  that  his  computer 
software  was  at  least  50  times  faster  than  IBM's  best  implementation  of 
Dantzig's  simplex  algorithm.  This  new  work  could  have  major  implications 
for  the  Casualty  Evacuation  Model  and  Kennington  spent  7  months,  full  time, 
investigating  this  new  algorithm.  This  investigation  is  reported  in 
Appendix  E.  The  group  found  that  while  the  new  method  clearly  worked,  it 
was  much  slower  than  their  existing  software  and  much  slower  than  the  claims 
of  Dr.  Karmarkar.  As  of  today,  no  group  has  been  able  to  substantiate 
the  claims  of  Dr.  Karmarkar.  Kennington  and  his  team  is  not,  at  this 
time,  investigating  the  projective  algorithm. 

While  Lt.  Col.  McLain  and  Capt.  Chmielewski  were  completing  the 
model  generator,  Professor  Kennington  and  a  former  student  (Mr.  David 
Allen)  worked  on  a  long-standing  problem  in  the  area  of  military 
communications  system  design.  The  problem  is  to  assign  frequencies  to 
nodes  to  minimize  both  co-channel  and  adjacent  channel  interference. 

They  designed  an  optimization  model  and  a  heuristic  solution  procedure 
which  virtually  solves  this  problem.  Problems  with  600  binary  decision 
variables  were  solved  in  only  10  seconds  on  a  moderately-sized  mainframe. 
This  work  is  reported  in  Appendix  F  and  this  paper  will  appear  in  Naval 
Research  Logistics  Quarterly. 

During  the  last  three  months,  Chmielewski  and  Kennington  have  been 
debugging  the  Casualty  Evacuation  Model  and  developing  the  appropriate 
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output  reports  from  Kennington's  solution  software.  The  work  so  far  has 
been  with  a  4-day  model.  The  goal  is  to  solve  a  full  90-day  model  during 
the  Spring  of  1986. 
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This  dissertation  presents  a  new  technique  for  solving  very  large 
scale  multicommodity  network  flow  problems.  The  method  obtains 
successively  better  lower  and  upper  bounds  on  the  optimal  objective 
function  value,  stopping  whenever  the  two  bounds  are  within  a 
prescribed  tolerance.  Lower  bounds  are  generated  by  partially  solving  a 
Lagrangian  dual  of  the  original  problem.  Upper  bounds  are  generated  by 
partially  solving  the  multicommodity  network  problem  itself,  using  a 
resource  directive  decomposition  scheme. 

One  great  advantage  to  our  approach  is  this:  in  both  the  lower  and 
upper  bound  routines,  the  problems  decompose  on  commodities.  As  a 
result,  only  a  single  commodity  minimum  cost  network  flow  optimizer  is 
needed.  Taking  advantage  of  this  natural  decomposition  also  allows  our 
technique  to  solve  substantially  larger  problems  than  other 
multicommodity  network  codes.  / 
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CHAPTER  I 
INTRODUCTION 

j 

This  dissertation  presents  a  new  technique  for  solving  very 
large  multicommodity  network  flow  problems. \  The  specific  application 
which  motivated  this  work  originated  with  the  United  States  Air  Force 
and  was  first  presented  to  us  by  Lt.  Col.  Dennis  McLain,  the  Assistant 
Director  of  Operations  Research  for  the  Military  Airlift  Command  at 
Scott  Air  Force  Base.  ^The  problem  is  an  extremely  large  casualty 
evacuation  model  to  be  used  by  the  Air  Force  in  Forming  a  plan  for  the 
evacuation  of  wartime  casualties.  This  plan  would  be  implemented  in 
case  of  a  European  military  conflict  involving  United  States  troops. 

SV-L 

Lt.  Col.  McLain  was  the  first  to  model  this  problem  as  a  multi- 
commodity  network  flow  problem  where  the  commodities  correspond  to  the 
various  types  of  wounds.  The  nodes  represent  such  entities  as  European 
bases  and  United  States  medical  facilities,  and  the  arcs  represent 
specific  aircraft  flights.  (A  complete  description  of  this  problem  is 
given  in  Section  1.3.)  This  problem  is  far  too  large  to  be  solved  by 
any  known  existing  computer  codes.  In  addition,  since  many  of  the 
data  are  only  rough  estimates  (the  number  of  casualties  of  various 
types  expected  at  given  locations),  an  exact  technique  is  not  cslled 
for.  Instead  a  technique  is  needed  to  discover  a  guaranteed  e -optimum 
for  any  given  c >0. 


This  is  precisely  what  our  technique  accomplishes.  It  generates 
successively  better  upper  and  lower  bounds  on  the  optimum,  stopping 
when  the  two  bounds  are  within  a  prescribed  tolerance.  We  exploit 
the  multicommodity  network  structure  in  both  the  lower  and  upper  bound 
routines  so  that  only  a  single  commodity  minimum  cost  network  flow 
optimizer  is  needed.  £VAC,  the  computer  code  which  implements  our 
technique,  has  been  used  to  solve  a  series  of  test  problems  in  less 
time  and  requiring  less  memory  than  MCNF,  a  specialized  multi- 
commodity  network  flow  problem  solver.  In  addition  EVAC  is  capable  of 
solving  very  large  problems  which  MCNF  is  unable  to  solve. 


1.1  Notation  and  Conventions 

The  notational  conventions  employed  throughout  this  work  are 

described  in  this  section.  Matrices  are  denoted  by  upper  case  Latin 

letters.  The  element  of  a  matrix,  A,  which  appears  in  the  ith  row  and 

jth  column  is  indicated  by  A^.  The  symbol  I  is  used  to  denote  an 

identity  matrix  with  dimension  appropriate  to  the  context.  Lower  case 

Latin  and  Greek  letters  are  used  to  denote  vectors.  The  symbol  0  is 

used  to  represent  a  vector  of  zeroes  with  dimension  appropriate  to 

the  context.  The  unit  vector,  whose  only  non-zero  component  is  a  one 

in  the  jth  position,  is  denoted  e..  Subscripts  are  used  to  indicate 

J 

individual  components  of  a  vector,  or  as  an  index  to  indicate  which  of 


3 


a  sequence  of  related  vectors  is  meant.  Superscripts  on  vectors  corre¬ 
spond  to  individual  commodities.  Note  that  vectors  are  considered  to 
be  row  vectors  or  column  vectors  as  appropriate  to  the  context;  that 
is,  no  special  notation  is  used  to  indicate  the  transpose  of  a  vector. 

The  inner  product  of  two  vectors,  x  and  y,  is  denoted  simply  by  xy. 

1  /2 

The  notation  ||x||  is  used  to  express  the  Euclidian  norm,  (xx) '  . 

Scalars  are  written  as  lower  case  Greek  or  Latin  letters. 

Euclidean  n-dimensional  space  is  denoted  Rn.  functions  are 
written  as  lower  case  Latin  letters,  and  functional  values  have  their 
arguments  in  parentheses.  For  example  g(y)  is  used  to  denote  the 
function  g  evaluated  at  the  point  y.  The  one  exception  to  this 
convention  is  the  projection  operation  described  in  Chapter  III.  In 
this  case  P[x]  is  used  to  express  the  projection  of  x  onto  the 
specified  region. 

Upper  case  Greek  letters  denote  sets,  with  the  exception  that 
8g(y)  is  used  to  denote  the  set  of  subgradients  of  a  function  g  at  a 
point  y  in  the  domain  of  g.  The  symbol  e  is  used  as  the  set  inclusion 
symbol  and  as  a  termination  tolerance. 

We  use  MAX{S}  to  denote  the  largest  element  of  a  set  S; 
similarly  MINIS)  indicates  the  smallest  element  of  S.  The  symbol  <*>  is 
used  for  infinity,  and  ■  denotes  the  end  of  a  proof.  All  other 
notation  is  standard. 
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1.2  Problem  Definition 

A  network  is  composed  of  two  entities:  nodes  and  arcs.  The  arcs 
may  be  viewed  as  undirectional  means  of  commodity  transport,  and  the 
nodes  may  be  thought  of  as  locations  or  terminals  connected  by  the 
arcs  and  served  by  whatever  physical  means  of  transport  are  associated 
with  the  arcs.  We  limit  our  consideration  to  networks  with  finite 
numbers  of  nodes  and  arcs.  For  a  given  network  we  denote  the  number 
of  nodes  by  m  and  the  number  of  arcs  by  n.  We  impose  an  ordering  on 
the  nodes  and  arcs  so  as  to  put  them  in  a  one-to-one  correspondence 
with  the  integers  {1,...,m}  and  {1,...,n},  respectively.  The  struc¬ 
ture  of  a  given  network  may  be  described,  then,  by  an  m  x  n  matrix 
called  a  node-arc  incidence  matrix.  Such  a  matrix,  A,  is  defined  in 
this  way: 

1+1 ,  if  arc  j  is  directed  away  from  node  i 
-1,  if  arc  j  is  directed  toward  node  i 

0,  otherwise. 

Additionally,  for  a  multicommodity  network,  we  are  concerned  with  more 
than  one  type  of  item  (commodity)  flowing  through  the  arcs.  We  order 
these  commodities  to  correspond  to  the  integers  {1 

We  define  the  following  quantities  to  be  used  in  the  formulation 
of  the  multicommodity  network  flow  problem: 

—  A  is  the  m  x  n  node-arc  incidence  matrix  corresponding  to  the 
underlying  network. 

—  x**  is  an  n  vector  of  decision  variables  for  k  =  1,...K.  Note 
k 

that  x.  represents  the  amount  of  flow  of  commodity  k  on  arc  j. 
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—  c  is  an  n  vector  of  unit  costs  for  k  s  1,...,K.  So 

c.  denotes  the  cost  for  one  unit  of  flow  of  commodity  k  on  arc 
J 

j. 

k  k 

—  r  is  an  m  vector  of  requirements  for  k  =  1,...,K,  so  that 

denotes  the  required  number  of  units  of  commodity  k  at  node  i.  If 

k 

<  0  then  node  i  is  said  to  be  a  demand  node  for  commodity  k 

with  demand  =  |rk|.  If  r*f  >  0  then  node  i  is  said  to  be  a 

1  i‘  l 

k  k 

supply  node  for  commodity  k  with  supply  =  r^.  And  if  r  :  0 

then  node  i  is  said  to  be  a  transshipment  node  for  commodity  k. 

—  u  is  an  n  vector  of  mutual  arc  capacities.  That  is,  the  total 

flow  of  all  the  commodities  combined  for  arc  j  cannot  exceed  u.. 

J 

—  v^  is  a  n  vector  of  arc  capacities  for  commodity  k  (k  = 

L- 

v  .  ,  then,  represents  an  upper  bound  on  the  flow  of  commodity  k 
J 

on  arc  j. 

We  sometimes  refer  to  the  entire  vector  of  decision  variables 
1  K 

(x  ,  •  • . , x  )  as  simply  x.  Similarly  we  use  c,  r,  and  v  to  denote  the 
entire  vector  of  costs,  requirements  and  upper  bounds,  respectively. 

Using  these  ideas,  we  may  formulate  the  multicommodity  network 
flow  problem  for  a  given  network  with  m  nodes,  n  arcs,  and  K  commodities 
as  follows: 

k  k 

Minimize  Sex 
k 

Subject  to  Ax^  =  r\  k  =  1,...,K  (MP) 

I  x^  u 
k 

k  *  1,...,K. 


n  .  k  .  k 
0  <  x  <  v  , 


1.3  The  Casualty  Evacuation  Model 


A  large  European  military  conflict  involving  U.S.  Armed  Forces 
could  result  in  more  casualties  than  could  be  effectively  handled  in 
European  medical  facilities.  To  alleviate  this  overcrowding,  the 
Department  of  Defense  plans  to  implement  the  following  evacuation 
policy: 

"During  the  first  30  days  of  a  conflict,  if  a  wounded 
soldier  cannot  be  returned  to  duty  within  15  days,  then 
he  will  be  evacuated  to  a  medical  facility  in  the  United 
States.  In  the  next  30  days  the  limit  on  treatment  time 
is  increased  to  30  days." 

Given  a  scenario  concerning  such  a  conflict  (i.e.  the  number  and  loca¬ 
tions  of  wounded  and  the  types  of  wounds) ,  this  evacuation  problem  may 
be  modelled  as  a  multicommodity  network  flow  problem.  Lt.  Col.  Dennis 
McLain  was  the  first  to  model  the  problem  in  this  way.  In  Lt.  Col. 
McLain's  model  the  nodes  correspond  to  9  European  recovery  bases  and 
95  United  States  locations.  The  arcs  correspond  to  aircraft  flights 
connecting  European  and  U.S.  facilities.  The  commodities  are  11 
different  patient  types. 

In  order  to  enforce  a  capacity  on  a  given  facility,  it  is 
necessary  to  duplicate  the  corresponding  node  using  the  capacity  as  an 
upper  bound  on  the  arc  between  the  duplicate  nodes.  For  example,  if 
node  A  represents  a  hospital  with  300  beds,  then  we  substitute  two 
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nodes,  A1  and  A2,  along  with  an  arc  whose  capacity  is  300.  Further, 
it  is  necessary  to  include  60  copies  of  the  entire  network,  one  for 
each  of  the  60  one-day  time  periods.  Additional  arcs  are  created  to 
link  each  time  period  to  the  next.  The  model  includes  a  dummy  "sink" 
node  for  each  time  period  and  one  "super  sink"  node,  along  with 
capacitated  arcs  to  allow  patients  who  have  recovered  to  exit  from  the 
system.  These  considerations  produce  a  very  large  model.  The 
dimensions  of  the  constraint  matrix  are  shown  below: 


A1 

*2 

• 

• 

• 

A.. 

11 

I 

I 

•  •  • 

I 

I 

12,561 


rows 


where  A^  =  ...  =  A^.  The  row  dimension  of  this  model  is  over  137,000, 
which  is  far  beyond  the  scope  of  any  known  existing  computer  code.  To 
put  these  figures  in  perspective,  we  note  that  Kennington  reports  that 
the  largest  models  he  has  solved  using  his  primal  partitioning  code, 
MCNF,  have  been  on  the  order  of  3000  rows  [2]. 

Our  plan  has  been  to  develop  a  specialized  solution  procedure 
which  would  solve  a  scaled-down  version  of  Lt.  Col.  McLain's  model.  We 
anticipate  aggregation  of  the  data,  possibly  using  some  of  the 
following  ideas: 
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(1)  Aggregation  of  the  time  periods.  Note  that  simply  using 
3-day  time  periods  instead  of  1-day  time  periods  reduces  the 
problem  size  to  around  46,000  rows. 

(2)  Aggregation  of  similar  patient  types. 

(3)  Aggregation  of  U.S.  medical  facilities  so  that  facilities 
which  are  located  within  a  given  number  of  miles  of  one 
another  are  treated  as  one  node. 

At  the  writing  of  this  dissertation  we  have  not  yet  received  any 
large  test  problems  from  the  Air  Force.  As  a  result,  we  are  unable  to 
report  on  the  problem  size  limitations  of  our  technique.  However,  in 
an  attempt  to  test  our  software  on  a  relatively  large  problem,  we 
solved  a  randomly  generated  test  problem  with  around  9,000  rows.  (See 
Chapter  4  for  the  details  of  this  problem.)  This  is  the  largest 
problem  we  have  attempted  so  far. 

1.4  Accomplishments  of  This  Investigation 

This  dissertation  proposes  a  new  technique  for  solving  extremely 
large  multicommodity  network  flow  problems.  Our  method  involves 
generating  upper  bounds  on  the  optimal  objective  value  by  partially 
solving  the  problem  using  a  resource-directive  decomposition  technique, 
and  generating  lower  bounds  on  the  optimal  objective  value  by  partially 
solving  a  Lagrangian  dual  of  the  problem.  Both  the  upper 
and  lower  bound  routines  exploit  the  network  structure  of  the  problem, 
decomposing  it  by  commodities  and  solving  the  resulting  pure  network 
problems.  In  the  limit  both  bounds  must  converge  to  the  optimal 
objective  function  value;  in  practice  we  stop  when  the  difference 
between  the  two  bounds  is  within  some  termination  tolerance. 


Whether  solving  for  lower  bounds  or  for  upper  bounds  a  sub- 
gradient  optimization  technique  is  used.  At  each  iteration  this 
procedure  requires  the  computation  of  a  subgradient,  the  selection  of 
a  step  size,  and  a  projection  operation.  In  Section  3.1,  we  obtain  a 
new  convergence  result  for  a  particular  class  of  subgradient  pro¬ 
cedures.  Then,  in  Section  3.2,  we  introduce  a  new  heuristic,  closely 
related  to  the  subgradient  optimization  procedure,  which  has  worked 
well  for  our  test  problems. 

Our  technique  has  been  tested  on  randomly  generated  test 
problems  and  on  one  problem  which  was  formulated  specifically  to 
represent  the  class  of  evacuation  planning  problems  for  which  the  code 
was  developed.  In  addition,  the  same  set  of  test  problems  was  solved 
by  MCNF  [51],  a  general  purpose  multicommodity  network  flow  problem 
solver  which  uses  a  primal  partitioning  scheme.  Computation  times  for 
both  codes  are  presented.  Our  code  used  an  average  of  68%  of  the  time 
needed  by  MCNF,  performing  significantly  better  on  the  problems  with 
fewer  commodities.  In  addition  our  code  required  on  the  order  of  1/K 
the  amount  of  main  memory  for  a  K-commodity  problem,  so  it  can  solve 
significantly  larger  problems  than  MCNF. 
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CHAPTER  II 

A  SURVEY  OF  RELATED  LITERATURE 


In  this  chapter  we  present  an  overview  of  the  existing  work  on 
which  this  dissertation  is  based.  Section  2.1  deals  with  the  work  that 
has  been  done  in  the  area  of  pure  network  models.  Then  in  Section  2.2 
we  address  the  broader  area  of  multicommodity  network  methods.  Since 
our  algorithm  involves  a  subgradient  optimization  technique,  both  in 
the  Lagrangian  dual  portion  and  in  the  resource-directive  decomposition 
routine,  we  provide  some  references  involving  subgradient  optimization 
in  Section  2.3 

2.1  Pure  Networks 

Network  problems  are  linear  programming  problems  with  node-arc 
incidence  matrices  as  their  constraint  matrices.  Within  this  class, 
known  formally  as  minimal  cost  network  flow  problems,  there  are  several 
variations  including  transportation  problems,  transshipment  problems, 
assignment  problems,  maximal  flow  problems,  snd  shortest  path  problems. 

Ideas  for  solution  of  network  problems  can  be  traced  at  least  as 
far  back  as  1939,  to  the  work  of  Professor  Leonio  Kantorovich  [41]. 
Kantorovich,  along  with  Professor  Tjalling  C.  Koopmans  received  the 
Nobel  Prize  in  Economic  Science  in  1975,  for  contributions  to  the 
theory  of  optimum  allocation  of  scarce  resources.  Koopmans  and  Reiter 
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[54]  and  Frank  L.  Hitchcock  [42],  working  independently,  were  the  first 
to  formulate  the  transportation  problem.  The  mid-fifties  saw  a  surge 
of  interest  and  work  in  the  areas  of  network  algorithms.  It  was  around 
this  time  at  Alex  Orden  [59]  generalized  the  transportation  model  to 
allow  transshipment  points.  Lester  Ford  and  Delbert  Fulkerson  [22] 

[20]  formulated  and  investigated  solution  techniques  for  the  maximal 
flow  problem  and  the  minimal  cost  network  flow  problem.  The  spe¬ 
cialized  algorithms  that  have  been  developed  for  solving  network 
problems  may  be  classified  into  two  groups:  primal-dual  techniques,  and 
specializations  of  the  primal  simplex  algorithm.  Primal-dual  methods 
for  solving  networks  began  with  Harold  Kuhn's  Hungarian  Algorithm  for 
the  assignment  problem  [55]  and  culminated  in  Fulkerson's  Out-of-Kilter 
Algorithm  [23].  Primal  simplex  based  techniques  originated  with  the 
work  of  Professor  George  Dantzig  [17]  and  continued  through  Ellis 
Johnson's  1965  paper  [47].  The  basis  for  Johnson's  work  can  be  traced 
to  the  work  of  Dantzig  [18]  and  Charnes  and  Cooper  [14], 

Since  that  time  much  work  has  been  done  in  the  area  of  solution 
techniques,  and  computational  advances  have  been  made  by  the  develop¬ 
ment  of  more  efficient  data  structures.  The  credit  for  much  of  this 
work  goes  to  Professors  Fred  Glover  and  Darwin  Klingman  and  their 
colleagues  at  the  University  of  Texas.  This  is  evidenced  by  such 
papers  as  Barr,  Glover  and  Klingman  [9]  [10],  Glover,  Hultz  and 
Klingman  [26]  [25],  Glover,  Karney  and  Klingman  [27],  Glover,  Karney, 
Klingman  and  Napier  [28],  Glover  and  Klingman  [29]  [31]  [30],  Glover, 
Klingman  and  Stutz  [32],  and  Karney  and  Klingman  [49].  Others  who  have 
contributed  to  the  research  include  Srinivasan  and  Thompson  [63]  [64], 
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Bradley,  Brown,  and  Graves  [13],  and  Mulvey  [57]  [58].  In  addition 
significant  work  has  been  performed  by  Professors  3eff  Kennington, 
Richard  Barr,  and  Richard  Helgason  of  Southern  Methodist  University  as 
seen  in  such  works  as  [3],  [41],  and  [52]. 

Today  network  algorithms  have  been  demonstrated  to  solve  linear 
network  problems  50  times  faster  than  general  linear  programming 
algorithms  [6].  Additionally  a  computer  implementation  of  such  a 
technique  may  require  only  half  the  memory  of  the  general  L.P.  package 
[6].  These  advances  are  due  to  the  efficient  data  structures  which 
have  been  developed  to  allow  a  basis  for  a  network  problem  to  be  stored 
as  a  rooted  spanning  tree  on  the  nodes  in  the  network.  Using  this  idea 
all  the  simplex  computations  such  as  pricing,  ratio  test,  and  updates, 
can  be  performed  via  labelling  algorithms  on  the  basis  tree.  This 
eliminates  the  need  to  store  the  basis  inverse  in  factored  form. 

2.2  Multicommodity  Networks 

Multicommodity  network  flow  problems  are  problems  in  which 
several  different  types  of  items  (commodities)  must  share  arcs  in  a 
capacitated  network.  Each  solution  technique  for  multicommodity 
network  models  can  be  classified  as  one  of  two  main  types  of 
algorithms:  partitioning  algorithms  and  decomposition  algorithms. 

2.2.1  Partitioning  Algorithms 

Partitioning  algorithms  are  specializations  of  the  simplex  method 
which  exploit  the  multicommodity  network  structure  by  partitioning  the 
basis  into  more  than  one  part.  In  one  part  advantage  is  taken  of  the 
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special  network  structure.  Those  who  have  studied  primal  partitioning 
algorithms  include  Kennington  [50],  Helgason  and  Kennington  [40],  Ali, 
Helgason,  Kennington,  and  Lall  [4],  Hartman  and  Lasdon  [36]  [35], 

Maier  [56],  and  Saigal  [61].  Ali  and  Kennington  [6],  in  their 
computational  research,  reported  solution  times  averaging  5  times 
faster  than  general  linear  programming  codes.  A  dual  partitioning 
method  was  proposed  by  Grigoriadis  and  White  [34].  A  primal-dual 
partitioning  scheme  was  developed  by  Oewell  [46].  In  addition  a 
factorization  technique  was  proposed  by  Graves  and  McBride  [33].  MCNF, 
the  multicommodity  network  code  with  which  we  compared  our  solution 
times,  is  a  primal  partitioning  program. 

2.2.2  Decomposition  Algorithms 

Decomposition  schemes  seek  to  solve  the  problem  by  decomposing  it 
into  several  smaller  subproblems,  each  of  which  takes  the  form  of  a 
pure  minimum  cost  network  flow  problem.  A  master  program  coordinates 
the  solution  process.  Decomposition  procedures  for  the  multicommodity 
network  flow  problem  fall  into  two  categories:  price-directive  schemes 
and  resource-directive  schemes. 

Price-directive  decomposition  is  based  on  the  well-known  research 
of  Dantzig  and  Wolfe  [19].  In  a  price-directive  approach,  the  K- 
commodity  problem  is  decomposed  into  K  single  commodity  problems.  The 
master  program  then  uses  the  simplex  method  while  the  subproblems  test 
for  optimality  and  select  candidates  to  enter  the  basis  of  the  master 
problem.  Ford  and  Fulkerson  [21]  were  the  first  to  develop  this 
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approach  for  solving  multicommodity  network  flow  problems.  Tomlin  [67] 
was  the  first  to  develop  a  computer  code  implementing  this  technique. 
Others  who  have  studied  price-directive  decomposition  schemes  are 
Darvis  [43],  Oarvis  and  Keith  [44],  Chen  and  DeWald  [15],  and  Oarvis 
and  Martinez  [45].  Price-directive  approaches  for  generalizations  of 
this  problem  have  been  proposed  by  Cremeans,  Smith,  and  Tyndall  [16], 
Swoveland  [65]  [66],  Weigel  and  Cremeans  [68],  and  Wollmer  [69]. 

Resource-directive  decomposition  schemes  decompose  the  problem  by 
commodities,  and  the  master  problem  systematically  distributes  the 
mutual  arc  capacities  among  the  commodities.  At  each  iteration  the 
optimal  solutions  to  the  single  commodity  subproblems  are  used  to 
compute  a  new  set  of  allocations.  Robacker  [60]  was  the  first  to 
suggest  this  approach  for  multicommodity  network  problems.  Research  on 
this  technique  has  been  presented  by  Swoveland  [65],  Assad  [8],  Ali, 
Helgason,  Kennington  and  Lall  [3],  and  Kennington  and  Shalaby  [53]. 

2.3  Subgradient  Optimization 

The  subgradient  optimization  technique  was  first  proposed  by  5hor 
[62]  in  1964.  Since  that  time  subgradient  algorithms  have  been  applied 
to  many  different  optimization  problems.  Held  and  Karp  [37]  and  Held, 
Wolfe  and  Crowder  [38]  made  use  of  the  approach  in  solving  the 
symmetric  travelling  salesman  problem.  Bazaraa  and  Goode  [11]  applied 
the  algorithm  to  the  asymmetric  travelling  salesman  problem.  Sub- 
gradient  methods  have  been  used  to  solve  the  assignment  problem  [38]. 
Glover,  Glover  and  Martinson  [24]  applied  a  subgradient  technique  to 
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solve  a  special  network  with  side  constraints,  and  Ali  and  Kennington 
[7]  made  use  of  it  in  research  involving  the  m-travelling  salesman 
problem. 
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CHAPTER  III 
THE  ALGORITHM 


Here  we  present  a  new  solution  technique  for  the  multicommodity 
network  flow  problem.  This  technique  involves  finding  successively 
better  upper  and  lower  bounds  on  the  optimal  objective  function  value. 
The  algorithm  terminates  whenever  the  two  bounds  are  within  a  prescribed 
tolerance  or  when  it  can  be  shown  that  the  current  solution  is  an  exact 
optimum. 

Lower  bounds  are  generated  by  partially  solving  a  Lagrangian  dual. 
At  each  iteration  a  Lagrangian  relaxation  of  the  original  problem  is 
solved;  since  these  relaxations  decompose  on  commodities,  only  a 
(single-commodity)  minimum  cost  network  flow  optimizer  is  needed.  A 
subgradient  direction  is  used  to  adjust  the  Lagrange  multipliers  for  the 
next  iteration. 

Upper  bounds  are  generated  using  a  modification  of  the  resource- 
directive  decomposition  technique  first  suggested  by  Robacker  [60].  We 
introduce  a  specialization  of  the  subgradient  direction  approach  which 
was  first  applied  to  this  class  of  problems  by  Held,  Wolfe,  and  Crowder 
[38]. 

With  minor  restrictions  on  the  step  sizes  we  show  that  both  the 
upper  and  lower  bounds  converge  to  the  optimal  objective  value  of  the 
original  multicommodity  network  flow  problem.  Hence  in  the  limit  the 
algorithm  will  converge  to  an  exact  optimum.  In  practice  we  seek  a 
near-optimum. 
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3.1  Subgradient  Optimization 

Let  us  first  consider  the  general  subgradient  algorithm  for 
optimization  of  convex  functions;  later  we  will  present  specializations 
of  the  technique  for  the  upper  and  lower  bound  problems.  Consider  the 
nonlinear  programming  problem 
Minimize  g(y) 

Subject  to  y  c  T 

where  g  is  a  real  valued  function  that  is  convex  over  the  compact, 
convex,  nonempty  set  F.  A  vector  n  is  called  a  subgradient  of  g  at  a 
point  x  if 

g(y)  -  g(x)  2  h(y  *  for  y  e  r  . 

Note  that  if  g  is  differentiable  at  x,  the  only  subgradient  at  x  is  the 
gradient.  We  denote  the  set  of  all  subgradients  of  g  at  x  by  3g(x). 

The  subgradient  algorithm  proceeds  in  this  manner:  Given  a  point 
x  in  r,  find  a  subgradient  of  g  at  x,  obtain  a  new  point  by  moving  a 
given  step  size  in  the  negative  subgradient  direction,  and  finally 
project  this  new  point  back  onto  r.  This  projection  operation  takes  a 
point  x  and  finds  the  point  in  r  that  is  "closest"  to  x  with  respect  to 
the  Euclidean  norm.  We  denote  the  projection  of  x  onto  r  by  P[x]. 

Using  this  notation  we  present  the  general  subgradient  optimization 
algorithm  for  minimizing  a  convex  function  g  [52]. 
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ALGORITHM  3.1  SUBGRADIENT  OPTIMIZATION  ALGORITHM 

Step  0  (Initialization) 

Let  be  any  element  of  r.  Select  a  set  of  step  sizes, 

S1,S2,S3’“*’  an<^  Se^ 

Step  1  (Find  Subgradient) 

Let  e  3g(y.).  If  n.  =  0  terminate  with  y^  optimal. 

Step  2  (Move  to  New  Point) 

Set  yA  Ptyj  -s^ni3.  Set  i i  +  1.  Return  to  step  1. 

Let  us  now  turn  our  attention  to  the  selection  of  step  sizes. 
Several  ideas  for  choosing  step  sizes  have  been  proposed.  These 
typically  involve  a  sequence  of  constants,  {a  ^ ,  A2 ,  A-j, . . .)  which  satisfy 
the  following  conditions: 

A  A  >  0,  for  all  i, 

lim  A^  =  0,  and  (3.1 ) 

i 

Z  A.  :  »  . 

1 

1 


The  subgradient  algorithm  can  be  shown  to  converge  when  any  of  the 
following  three  formulae  are  used  for  determining  step  sizes  [52]: 


(ii)  sA  =  Aa/| | nA ! | 2  *  (3.2) 

(iii)  sA  =  Ai[g(yi)  -  g*3/||ni||2 


where  g*  denotes  the  optimal  objective  value. 
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Propositions  3.1,  3.2,  and  3.3  may  be  found  in  Kennington  and 
Helgason  [52],  and  are  given  here  as  necessary  preliminary  results. 
Proposition  3.1  [52] 

Let  yeT,  and  let  xcRn.  Then  (x-P[x])(y-P[x])  £0. 

Proof 

Choose  a  so  that  0<a<1 .  Since  r  is  convex,  ay+(1-a)P[x] er.  By 
the  definition  of  P[x],  | |x-P[x] |  )£| |x-(ay+(1-a)P[x]) | | .  Thus 
| |x-P[x]||2  <  | |x-(ay+(1-a)P[x]) | |2 
=  | |x-P[x]-a(y-P[x]) | |2 

=  | |x-P[x] | j2+a2| |y-P[x] | |2  -2a(x-P[x])(y-P[x]) . 

Then  (x-P[x])(y-P[x])  £  | | y-P[x] | | a/2 .  And,  since  a  can  be  taken 
arbitrarily  close  to  0, 

(x-P[x] )(y-P[x] )  £  0  .  ■ 

Proposition  3.2  [52] 

Let  x,  y  e  Rn.  Then  | |P[x]-P[y] | |  £  ||x-y||. 

Proof 

Case  1:  Supose  P[x]  =  P[y].  Then 
| |P[x]-P[y] j |  =  0  £  | | x— y | | . 

Case  2:  Suppose  P[x]  i  P[y].  Then  since  P[x]cF, 
and  P[y]eT,  from  Proposition  3.1  we  have  that 
(x-P[x])(P[y]-P[x])  £0 
and 

(y-P[y])(P[x]-P[y])  £  0. 

We  may  rewrite  the  above  inequalities  as 
x(P[y]-P[x])-P[x]P[y]+j |P[x]||2  £  0 


and 


y(P[x]-P[y])-P[y]P[x]-*-j  |P[y]|  l2  <  0. 

Adding  these  inequalities,  we  obtain 
(x-y)(P[y]-P[x])+| |P[y3-P[x] | |2  <  0. 

Then  from  the  Cauchy-Schwartz  inequality, 
-(x-y)(P[y]-P[x])  <  | | x-y | |  | |P[ y]-P[ x] | | . 
Thus 

|  |P[y]-P[x]  1 12  <  | | x-y | |  I |P[y]-P[x]  |  I . 


And  since  P[x]  i  P[y], 

||P[x]-P[y]||  <  ||x-y||.  . 

Proposition  3.3  [52] 

If  i  0,  then,  for  any  yeT, 

Mvu1-y||2  <  llyi-vll2  -  •12llnill2-*B1'.i<y-y1). 
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Proof 

Let  i  be  any  iteration  of  the  subgradient  algorithm.  Suppose 
/  0.  Let  yeT.  Then,  by  Proposition  3.2, 

I  Ipt>,i-sirii3-pty3 }  |2  <  I  |yi-sini-y|  |2 


=  I  Wi-y  1 12  +  s^llnjl2  +  2sini(y-yi) . 

Since  P[y]=y  and  P[y^-s.n^J  =  y^,  we  *"iave  that 


l|yM-yll2  <  llyi-yll2  -  si2ll"ill2  *  zvi'y-yi’  •  ■ 

Our  main  convergence  result  is  for  the  particular  step  size 
scheme : 

Si  =  Xit9<yi>  -  g]/||„il|2 


t 
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where  g  is  a  lower  bound  for  the  optimal  objective  and  where  we  are  at 
liberty  to  select  bounds  a  and  6  tor  the  such  that  for  each  i,  0  < 

a  <  Xi  <  6  <  2. 

Proposition  3.4 

Let  (i)  g  be  a  known  lower  bound  for  the  optimal  objective,  g*, 
with  g*>g; 

(ii)  {\.}  be  any  infinite  sequence  such  that 
for  all  i,  0<aO^<_8<2;  and 
(iii)  si  =  Aiig(yi)-g]/| | n± | |2- 

If  there  is  a  constant  C  such  that  for  all  i,  Ifn^ll  £  C,  and  if  y  >  0 
is  given,  then  there  is  some  n  such  that  g( yp )  £  g*+[ e/(2-B)](g*-g)+ 

y 

Proof 

Let  y>0  be  given.  Let  (i),  (ii),  and  (iii)  hold.  Let  y*  be  an 
optimal  point,  and  for  all  i,  |  | n ^ | j  £  C.  Suppose,  contrary  to  the 
desired  result,  that  for  all  n,  g(yn)>g*+[  6/(2-B)](g*-g)+'y.  Then,  by 
Proposition  3.3, 

lly^y-ll2  i  l|yi-y||J-x12[»(yi)-5]2/|U1ll2 

*  21-i it gC  Vj >-Q]/ 1  Irij  I  ihr^ty'-Yj) 

<  l|yi-y||2*iiZ[g(y1)-g]2/||nil|Z 

*  2xi{[g<yi)-gl/|  IgJ  |2  )Ig*-g(yi'], 
since  nic9g(yi). 

Since  B£Xi>0,then  BX^  X^.  So, 


! Iyi+ry*' !2  -  Myi-y*ll2+6>'it9(yi)-g]2/||ni||2 
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+  2Xi{[g(yi)-g]/| | n± | |2 }[g*-g(yi) ] 

=  I  |yi-y*M2+(2-B)xi{tg(y1)-g]/||ni|  I2} 

[ ( g*-g( y p )+( b/(2- b) ) ( g*-g) ] . 

Since  gC y ± )  >  g*+( B/(2-s) ) ( g*-g)+Y ,  then  -Y  >  g*-g(yi)  + 
(B/(2-e))(g*-g) •  so, 

I  lyit-ry*l  I2  <  1 1  y  i~y*  1  i2-(2-B)xi[g(yi)-g]Y/|lnil  I2. 

Since  g*_<g(yA) ,  a£X^  ,  and  ||n^|(  £  C,  then 

l|yi+1-y*ll2  <  l|yi-y*||2-[(2-6)a(g*-g)Y]/c2.  (3.3) 

We  can  choose  an  integer  N  so  large  that 
C2||y1-y*||2/(2-6)a(g*-g)Y  <  N. 

Thus,  since  2-B>0  and  g*-g>0, 

N(2-B)a(g*-g)Y/C2  >  ||yi-y*||2. 

Adding  together  the  inequalities  obtained  from  (3.3)  by  letting  i  take 
on  all  values  from  1  to  N,  we  obtain 

l|yN+ry*ll2  <  ||yry*||2-N(2-B)a(g*-g)Y/c2  <  o, 


a  contradiction.  ■ 
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It  is  shown  in  [39]  that  when  r  is  compact,  g  is  continuous  on  some  open 
set  containing  r,  and  3g( y )  i  4,  for  all  yer,  there  exists  a 
constant  C  such  that  ||n||_£C  for  all  ycT*  and  ne39(y)»  so  that  the 
boundedness  condition  on  the  subgradients  in  Proposition  3.4  is  easily 
met. 


3.2  Generating  Lower  Bounds 

In  this  section  we  present  a  technique  for  generating  lower  bounds 
for  the  multicommodity  network  flow  problem.  This  technique  involves 
partially  solving  the  Lagrangian  dual  problem  using  a  subgradient 
technique  to  update  the  Lagrange  multipliers  at  each  iteration. 

Recall  that  the  multicommodity  network  flow  problem,  MP,  may  be 
stated  as  follows: 

k  k 

Minimize  I  c  x 
k 

Subject  to  Axk  =  rk,  k  =  1,...,K  (MP) 

I  xk  <  u 

k 

0  <  xk  i  vk,  k*1,...,K 


where 

A  is  an  m  x  n  node-arc  incidence  matrix, 

|( 

c  is  an  n  vector  of  unit  costs  for  k  =  1,...,K, 

k 

r  is  an  m  vector  of  node  requirements  for  k  =  1,...,K, 
u  is  an  n  vector  of  mutual  arc  capacities, 

k 

v  is  an  n  vector  of  individual  commodity  bounds  for  k=1,...,K, 
x  is  an  n  vector  of  decision  variables  for  k  = 

, . . . ,  K , 


1 
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and  K  is  the  number  of  commodities. 

Consider  a  Lagrangian  dual  problem  for  MP,  denoted  by  DP: 

MAX  h( a) 

A_>0 

h( A)  =  MIN[  I  ckxk  +  A(I  xk  -  u):  (DP) 

k  k 

Axk  =  rk  (k  =  1 ,...,«)?  0  <  xk  <_  vk  (k  =  1 . K)] 

where  A  is  an  n  vector  of  Lagrange  multipliers. 

First  we  show  that  any  feasible  solution  for  DP  is  a  lower  bound 
for  MP. 

Proposition  3.5  [12] 

-  -1  -2  -K 

Let  x  =  (x  ,x  , ...,x  )  be  a  feasible  solution  for  MP.  Let 
A  be  a  feasible  solution  for  DP.  Then  h(X)  £  cx. 

Proof 

Since  h(A)  is  a  minimum,  and  since  x  is  feasible  for  MP,  h(X)< 

k  __  — ___ 

Zc  x  +  A(lx  -  u) .  Further  since  A  is  feasible 

k  k  -  -k 

for  DP  and  x  is  feasible  for  MP,  then  A(Zx  -  u)  £  0. 

_  k 

Hence  h(A)  <_  cx.  ■ 

In  addition  to  this  result,  Bazaraa  and  Shetty  [12]  have  proved 
that  if  MP  has  an  optimal  solution,  then  DP  has  an  optimal  solution,  and 
that  their  optimal  objective  function  values  are  equal.  As  a  result,  we 
see  that  we  may  indeed  solve  (or  partially  solve)  DP  in  order  to  obtain 
a  lower  bound  for  MP. 

In  order  to  justify  using  a  subgradient  optimization  technique  for 
solving  DP,  we  must  show  that  the  objective  function  is  concave  and 
develop  an  expression  for  a  subgradient. 


2b 


Proposition  3-6 

The  real  valued  function  h  is  concave  over  A  =  (X:X  cRn;  X  >  0} 


Proof 

Let  X1  >  0.  Let  X2  >0.  Let  Q  <  a  <  1  •  Then 

hUX*  +  (1-a)  X2 )  =  MINI  I  ckxk  +  (aX1  ♦  (l-a)X2)(Ex  -u): 

Axk  =  rk(k=T,...,K);0  iKk  £vk  (k=1,...,K)] 

-  MIN[aIckxk  +  aX1  (Zxk-u)+(1  -o)Eckxk  +  (1- a)  X2(Exk-u) : 
k  k  k 


Axk  = 


:(k=T,...,K);  0  <  x  <  ' 


‘(k=1 


,K)J 


>  aM!N[Ickxk  +  X^lx  -u): 

k  k 


Axk  s  rk(k=1  .  ,K);0  <  xk  <_  vk(k=1  , . . .  ,K)  ] 

+  (1-a)  MIN  [ickxk  +  X2(Exk-u)  : 
k  k 

Axk  =  rk(k=1  ,...,K);0  (k=1  , . . . ,«)] 

=  ahU1)  +  (l-a)h(X2).  Hence  h  is  concave  over  A. 


Proposition  3.7 

Let  \  >  0*  Let  x  represent  an  optimal  value  of  x  correspon  g 
to  hOO.  Then  d  =  j*k-ti  is  a  subgradient  of  h  at  X. 

Proof 

Let  \  be  any  other  point  in  A  -ith  corresponding  optimal  decision 

A 

variable  values  x.  Then 

h(x)  =  I  ckxk  +  XU  xk-u) 
k  k 
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k-k  *  _k  * 

<_  Z  c  x  +  X(Z  x  -u)  (since  x  is  optimal) 

k  k 

=  I  ck5k  +  X(  Z  xk-u)  +  (I  Z  xk-XZ  xk)  +  (lu-^u) 

k  k  k  k 

=  Z  ckxk  +  1(Z  xk-u)  +  (Z  xk-u)(X-3.) 
k  k  k 

=  h(X)  +  d(  X- X) . 

Therefore  d  is  a  subgradient  of  h  at  X.  ■ 

We  now  present  our  algorithm  for  computing  lower  bounds  for  MP. 
Note  that  it  is  a  specialization  of  the  subgradient  optimization 
algorithm  for  this  problem,  and  its  convergence  follows  as  a  maximiza¬ 
tion  analog  of  Proposition  3.4. 

ALGORITHM  3.2  LOWER  BOUND  ALGORITHM 


Step  0  (Initialization) 

Let  UB  be  any  upper  bound  on  the  solution  to  MP.  Set  i  0; 

1  K 

0;  aQ  +  2.  Compute  yQ  +•  h(XQ)  and  let  xQ  =  (x0,...,xQ)  be  the 

corresponding  optimal  values  of  the  decision  variables. 

Step  1  (Find  Subgradient) 

1/ 

Set  n.  +  Z  x.-u.  If  n-  =  0,  stop  with  optimal, 
k 

Step  2  (Move  to  New  Point) 

2  th 

Set  +•  cu  (UB-y^)/| | r^| |  •  Compute  the  j  component  of 


X.  , 

i+l 


as: 


(xi+1}j 


MAX  {(xi+sini)j  ,  0} 
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Compute  •*-  h(  )  and  let  x^  be  the  corresponding  optimal  values 

of  the  decision  variables.  Set  a-/2.  Set  i  i+1 .  Return  to 

a+1  1 

step  1 . 

3.3  Generating  Upper  Bounds 

Here  we  describe  a  procedure  for  generating  upper  bounds  for  the 
multicommodity  network  flow  problem.  This  procedure  is  a  specialization 
of  the  resource-directive  decomposition  (ROD)  algorithm  using  a  sub¬ 
gradient  direction.  First  we  describe  the  general  RDD  procedure;  then 
we  present  our  specialization. 

The  RDD  technique  produces  a  sequence  of  feasible  solutions  by 
distributing  the  mutual  arc  capacity  among  commodities  in  such  a  way 
that  the  solutions  to  the  K  individual  subproblems  provide  a  solution  to 
the  composite  problem.  At  each  iteration  an  allocation  is  made  and  the 
resulting  K  (single  commodity)  minimum  cost  network  flow  problems  are 
solved.  If  the  solution  meets  an  optimality  criterion  then  the 
procedure  terminates;  otherwise,  a  new  allocation  is  made,  and  the 
process  is  repeated. 

After  introducing  artificial  variables,  (a  ),  MP  becomes: 

Minimize  E  ckxk  +  M  I  1ak 
k  k 

Subject  to  Ax'*"  +  ak  =  rk  (k  r 
_  k 

Ex  <_  u 
k 

0  <_  xk  <_  vk  (k  :  1,...,K) 
ak  >_  0  (k  =  1 , . . .  ,K) 
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where  M  is  a  very  large  positive  number  and  1  is  an  m  vector  of  all 
ones. 

Let  us  restate  the  problem  as: 

1  K 

Minimize  z(y  , . . . ,y  ) 

Subject  to  zCy1 , . . . ,yK)  =  I  zk(yk)  (RP) 

k 

„  k 

Zy  =  u 
k 

0  <yk<  vk  (k  =  1 , . . .  ,K) 

1/  1/  k  L-  Lf  L/  1/ 

where  z  (y  )  =  MIN  {c  x  -*-M_1_a  :  Ax+a  =  r  ;  Ckx^y  ;  a>0}  for  k  = 

k  k 

We  shall  refer  to  this  formulation  as  RP.  Note  that  z  (y  )  = 

...  v  r  k  k  k  k  k.  k  k  k  ...  kni  .  ,  ... 

MAXir  p  -y  v  :p  A-v  <_c  ;p  <^M_1_;v  >0),  by  duality  theory. 

In  order  to  justify  using  a  subgradient  optimization  technique  we 
1  K 

must  show  that  z(y  ,...,y  )  is  a  convex  function  and  develop  an 
expression  for  a  subgradient. 

Proposition  3.8  [52] 

The  real  valued  function  z  is  convex  over 
Y  =  i y  , . . . ,  y  : y  2. 0 ; . . . ; y  >0}. 

Proof 

Let  (y\...,y*^)eY  and  (y\ . . .  ,yK )cY.  Select  a  so  that 
0<cx<1-  Then 

z[ay  *K1-ct  )y  ,...,ay  +  (1-a'y  J 

kr  -k  ,‘k-, 

=  r  z  [ay  +(l-a)y  ] 

k 

=  I  MAX  rkuk-  [ayk  +  (1-a)yk]vk: 
k 

k.  k  k  k  ...  k  n. 
p  A-v  <C  ;  p  <  Ml  ;v  >0} 
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iiav  <  f  ^  ^  /  . f  xrkk^kk-* 

=  I  MAXialr  y  -y  v  ]  +  (1-a)[r  y  -y  v  J: 

k 

k.  k  k  k  ...  k„ 
y  A-v  <c  ;  y  <M_1_;  v^_0} 

£  a I  MAX  { rkyk-ykvk: 

k 

k.  k  k  k  ...  k  n  , 
y  A -v  <c  ;  y  <M1  ;  v  >0  } 

+  ( 1  — a )  I  MAX{  rkyk-ykyk: 

k 

k.  k  k  k  ...  k 

y  A-v  _<c  ;  y  v  >0} 

1  -Kx  \  r\  *Kv 

=  az(y  ,...,y  )  +  (l-a)z(y  ,...,y  ). 

Therefore  z  is  convex  over  Y.  ■ 


Proposition  3.9  [52] 

-  -1  -K  -k  -k 

Let  y  =  (y  ,...,y  )cY  be  any  allocation  and  let  ( y  ,v  ) 

denote  the  corresponding  optimal  solution  to  z  (y  )  for  k  = 

-1  -K  - 

Then  n  =  (-v  ,  )  is  a  subgradient  of  z  at  y. 

Proof 

IK  k  k 

Let  y  =  (y  ,...,y  )cY  be  any  allocation  and  let  (y  ,v  )  denote  the 

k  k 

corresponding  optimal  solution  to  z  (y  )  for  k  =  Then: 

/  1  Kx  ,-1  -K\  r  ,  k  k  k  ks  _  ,  k-k  -k-kx 

z(y  ,...,y  )-z( y  ,...,y  )  =  l  (r  y  -y  v  )  -  l  (r  y  -y  v  ) 

k  k 

.  ,  k-k  k-kx  _  ,  k-k  -k-kx 

>_  l.  (ry-yv  )  -  I  (ry-yv  ) 


=  I  (_:k)(yk_yk)  . 
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Hence  is  a  subgradient  of  2  at  y.  ■ 

Recall  that  the  subgradient  optimization  algorithm  requires  a 

technique  for  projecting  a  point  onto  the  feasible  region.  We  now 

explore  the  projection  operation  for  this  problem. 

Let  us  denote  the  feasible  region  for  RP  by  0 .  That  is, 

ft  =  {(y1 , . . .  ,yK) :  £yk  =  u;0  5  yk  ±  vk(k  =1 , . . . ,K) }.  Given  an 

-1  -K 

arbitrary  allocation,  (y  ,...,y  ),  to  project  it  onto  ft,  we  solve  : 

MIN{ | | (y1 ,. .. ,yK)-(y1 ,. . . ,yK)  |  j :  y  eft} 

=  MIN  {(ll(yk-yk)2)1/2  :  y  eft}, 
kj  J  J 

Or,  equivalently,  we  can  solve: 

MIN{n;(yk-yk)2  :  y  efi}. 
kj  J  J 


Note  that  this  problem  decomposes  on  j.  Hence,  for  each  arc  j,  we 
solve: 


MIN(l(yk-yk)2  :  I  yk  =  u  ;  0  <  yk  <  vk  ;  (k=1 , . . . ,K)  }. 


We  will  denote  the  above  projection  problem  by  P.  The  following 
algorithm  [52]  is  used  to  solve  P  for  any  arc,  j. 


ALGORITHM  3.3  PROJECTION  ALGORITHM 


Step  0  (Initialization) 
u 

If  u.  >  I  v.  or  u.  <  0,  terminate  with  no  feasible 

J  k  J  J 


solution.  Otherwise  set  1  1 ;  r-<-2K;  L-*-I  v  R-H3.  Compute 

k  J 
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-k  -k  k 

the  breakpoints,  b.  (i=1 , . . . ,2K),  as  y.  and  y.-v.  (k=1,...,K). 

*  J  J  J 

Order  the  breakpoints  so  that  b^  <_ 


If  r-1  =1  go  to  step  4;  otherwise,  set  m<-[(l+r)/2]j  where  [K]^ 

is  the  greatest  integer  <  k. 

Step  2  (Computer  New  Value) 

Set  C  +  j  MAX(MIN[yk-y  ,  vk],  0} 
k  J  m  J 

Step  3  (Update) 

If  C=c  then  set  A-<-y*  and  go  to  step  5.  If  C>c  then 

m 

set  l-«-m:  L-*-C;  and  go  to  step  1.  If  C<c  then  set  r-^m;  R-*£;  and  go  to 
step  1 . 

Step  4  (Interpolate) 

Set  A*^b.+[(b  -b.)(c-l)]/(R-l). 
i  r  1 

Step  3 

Compute  the  feasible  (projected)  allocation,  y.,  for 

J 

k=1,...,K  in  this  way: 


Ik  .c  -k  k 

v j  ’  if  X*  1  Yj-Vj 

-k  , _  . c  -k  k,  _  -k 

y  -A*  ,  if  y.-v.<A*<y- 

J  J  J  ”  J 

0  ,  if  A*  > 

J 


Terminate  with  the  feasible  allocation  for  arc  j, 


An  upper  bound  algorithm  using  the  subgradient  procedure  is  now 


presented.  Its  convergence  is  a  direct  result  of  Proposition  3.4. 


32 


ALGORITHM  3.4  UPPER  BOUND  ALGORITHM 

Step  0  (Initialization) 

Let  LB  be  any  lower  bound  on  the  solution  to  MP.  Choose  a  set  of 

IK  k 

initial  allocations,  y^  =  (yg,...,yg)  by  setting  y^  P[  (1/K)(u)  ] 
for  k  =  Set  Xg  •*-  2;  i  ■*-  0;  UB-«o. 

Step  1  (Find  Subgradient) 

Let  (pk,vk)  solve  zk(yk)  for  k  =  Let 

ni"'(”vi ’  *  *  *  »"v^)  •  Set  UB-t-zk(yk).  If  =  0,  then  terminate  with 
z(y^)  optimal. 

Step  2  (Move  to  New  Point) 

Compute  si  •<-  Xitz(yi)-LB]/j  |ni|  |2.  Set  y.+1  «. 

P^yi_sini^‘  5et  Xi+1^i//2;  *  Return  to  step  1  • 

We  now  introduce  a  heuristic  modification  of  the  upper  bound 

algorithm,  which  has  produced  better  results  on  our  test  problems. 

IK  IK 

Recall  that  n  =  (-v  ,...,-v  )  is  a  subgradient  of  z  at  (y  ,...,y  ). 

Then  for  each  arc  j,  the  vector 

n(j)  =  (h.e.,  g  .  e  . ,. ..  ,n/.  „\  .  e.) 

J  J  j’  m-j  j’  ’  (k-1  )n+j  y 

serves  to  isolate  the  components  of  n  associated  with  the  commodities 

flowing  on  arc  j.  For  each  such  arc  j  we  compute  an  individual  step 

size  at  iteration  i  as 

si(j)  =  Aitz(y1.,...,y*)-z*]/|  |n  ±(  j)  1 12 


where  z*  is  approximated  by  LB. 
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Using  this  idea  we  now  present  our  heuristic  upper  bound 
algorithm. 


ALGORITHM  3.3  HEURISTIC  UPPER  BOUND  ALGORITHM 

Step  0  (Initialization) 

Let  LB  be  any  lower  bound  on  the  solution  to  MP.  Choose  a  set  of 

IK  k 

initial  allocations,  y^  =  (y^,...^)  by  setting  y^  *■  P[(1/K)(u)] 

for  k  =  Set  Ag«-2;  i-*-0;  UB-«°. 

Step  1  (Find  Subgradient) 

Let  (pk,vk)  solve  zk(yk)  for  k  =  Let  = 

(-vl , . . .  ,-vf) .  Set  UB-*-Lzk(yk) .  If  rK=0,  then  terminate  with 
li  k  i 

z(y,^)  optimal. 

Step  2  (Move  to  New  Point) 

IK  2 

Compute  si(j)-t-ii[z(yi,...  ,yi)-LB]/|  |ni(j)  1 1  for  each  arc  j. 

Set  S  *-  diag(si(1 ) , . . .  jS^n) ) .  Set 


Set  (yj+1  Set  xi+i^xi/2;  i"Hi+1* 

Go  to  step  1 . 


In  this  section  we  present  the  composite  algorithm  for  solving  MP. 


This  procedure  involves  partially  solving  DP  for  successively  better 
lower  bounds  and  partially  solving  RP  for  successively  better  upper 
bounds  on  the  optimal  objective  function  value.  The  algorithm 
terminates  whenever  (a)  the  solution  to  DP  can  be  shown  to  be  an  exact 
optimum;  (b)  the  solution  to  RP  can  be  shown  to  be  an  exact  optimum;  or 
(c)  the  greatest  lower  bound  and  the  least  upper  bound  generated  are 
within  a  prescribed  tolerance,  e  .  In  case  (c),  the  best  solution  to  RP 
is  presented  as  a  guaranteed  e-optimal  solution. 


ALGORITHM  3.6  COMPLETE  ALGORITHM 

Step  0  (Initialization) 

Let  extermination  tolerance  (0<e<1);  NDLB-xnumber  of  lower  bound 
iterations  to  perform  on  each  pass;  NDUB+number  of  upper  bound 
iterations  to  perform  on  each  pass;  LB-*--®;  UB-*-®. 

Step  1  (Lower  Bound) 

Perform  NOLB  iterations  of  the  lower  bound  algorithm  (Algorithm 
3.2).  Let  LB  denote  the  best  lower  bound  attained  so  far.  If  Algorithm 
3.2  terminates  in  step  1  with  an  exact  optimum,  terminate  with  that 
solution  optimal  for  MP. 

Step  2  (Upper  Bound) 

Perform  NOUB  iterations  of  an  upper  bound  algorithm  (Algorithm  3.4 
or  3.5).  Let  UB  denote  the  best  upper  bound  attained  so  far.  If 
Algorithm  3.4  terminates  in  step  1  with  an  exact  optimum,  terminate  with 
that  solution  optimal  for  MP. 
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Step  3  (Check  for  Termination) 

If  e(UB)<LB  then  terminate  with  UB  a  guaranteed  e-optimum; 
otherwise,  go  to  step  1. 

In  this  algorithm  the  best  solutions  for  the  lower  bound  and  upper 
bound  problems  at  each  pass  are  retained  and  used  as  starting  solutions 
for  the  respective  problems  on  the  next  pass.  The  details  of  our 
implementation  are  presented  in  Chapter  4. 
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CHAPTER  IV 

COMPUTATIONAL  EXPERIMENTATION 


This  chapter  provides  descriptions  of  our  computer  implementation 
of  Algorithm  3.6  and  of  the  test  problems  used.  Our  code,  EVAC,  uses 
MODFLO  [l]  to  solve  the  single  commodity  minimum  cost  network  flow 
subproblems  which  arise  in  Algorithm  3.2  and  in  Algorithm  3.5.  MODFLO 
is  a  set  of  routines  which  may  be  used  to  solve  a  network  flow  problem  or 
to  reoptimize  a  previously  solved  problem  after  changes  are  made  in  some 
of  the  data.  MODFLO,  which  is  based  on  NETFLO  [52],  allows  the  user  to 
change  bounds,  costs,  and/or  requirements  and  then  reoptimize  from  a 
basis  which  was  optimal  for  the  original  problem. 

We  tested  EVAC  on  22  randomly  generated  multicommodity  network 
flow  problems  and  on  one  test  problem  which  was  specially  structured  to 
be  solved  by  EVAC.  The  test  problems  ranged  in  size  from  22  to  754 
nodes  and  from  53  to  1 ,102  arcs  with  from  0  to  599  linking  constraints 
and  from  3  to  20  commodities.  The  equivalent  LP  sizes  are  between  232 
and  8,904  rows  and  between  470  and  12,111  columns.  The  22  randomly 
generated  problems  were  created  using  MNETGN  [5],  a  multicommodity 
network  problem  generator.  The  problems  were  solved  by  EVAC  and  by  MCNF 
[51],  a  multicommodity  network  flow  code  which  uses  a  primal  parti¬ 
tioning  algorithm.  Solution  times  are  compared  and  conclusions  are 
drawn  concerning  the  relative  effectiveness  of  the  techniques. 
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4.1  Description  of  the  Computer  Programs 

In  this  section  we  present  a  description  of  MCNF  and  EVAC,  the  two 
computer  codes  used  in  our  experimentation.  Both  programs  are  written 
in  standard  FORTRAN  and  have  been  tailored  to  neither  our  equipment  nor 
our  FORTRAN  compiler. 

4.1.1  MCNF 

MCNF  was  developed  by  Oeff  Kennington  at  Southern  Methodist 
University,  Dallas,  TX.  It  is  an  incore  multicommodity  network  flow 
problem  solver  which  uses  the  modification  of  the  revised  simplex  method 
known  as  the  primal  partitioning  algorithm  [36].  In  this  algorithm  the 
basis  inverse  is  maintained  as  a  set  of  rooted  spanning  trees  (one  for 
each  commodity)  and  a  working  basis  inverse  is  maintained  in  product 
form.  The  working  basis  inverse  has  dimension  equal  to  the  number  of 
binding  linking  constraints  corresponding  to  the  current  basis.  The 
initial  basis  is  created  using  a  multicommodity  variation  of  the  routine 
used  in  NETFLO.  A  partial  pricing  scheme  is  used;  the  pricing  tolerance 
is  1.E-6  and  the  pivot  tolerance  is  1.E-8. 

4.1.2  EVAC 

EVAC  is  our  implementation  of  Algorithm  3.6  for  solving  the 
multicommodity  network  flow  problem.  Note  that  Algorithm  3.6  alternates 
between  generating  lower  bounds  using  Algorithm  3.2  and  generating  upper 
bounds  using  Algorithm  3.5.  Since  both  the  lower  bound  problem  (DP)  and 
the  upper  bound  problem  (RP)  decompose  on  commodities,  EVAC  maintains 
only  the  information  concerning  the  current  commodity  in  main  memory. 

The  problem  data  and  most  recent  bases  for  all  the  other  commodities  are 
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kept  on  peripheral  storage.  At  the  user's  option  EVAC  stores  in  main 

1  k 

memory  as  much  of  the  current  set  of  allocations,  (y^,...,y^)  and 

1  k 

current  dual  variables  (-v^, . . .  ,-v^)  as  desired.  All  our  test 
problems  (with  the  exception  of  Problem  23)  were  solved  with  all  the 
allocations  and  dual  variables  in  core. 

Both  the  lower  bound  routine  and  the  upper  bound  routine  use 
MODfLO  as  the  optimizer  for  the  single  commodity  subproblems.  MODFLO 
uses  the  same  partial  pricing  scheme  as  NETFLO  and  drives  the  flow  on 
artificial  arcs  to  zero  using  the  Big-M  method.  The  Big-M  value  that 
was  used  for  our  test  problems,  except  as  noted  in  Table  4.1,  was  7 
times  the  largest  unit  cost  in  the  given  problem.  At  subsequent 
iterations,  initial  bases  for  each  commodity  are  just  the  optimal  bases 
for  the  previous  set  of  Lagrange  multipliers.  A  basis  for  the  upper- 
bound  problem  is  generated  by  constructing  a  feasible  basis  from  the 
previous  optimal  basis  using  the  rules  described  in  [13- 

In  practice  we  did  not  update  the  multipliers  for  the  step  sizes 
( cu  in  Algorithm  3.2  and  in  Algorithm  3.5)  at  every  iteration,  but 
only  when  the  improvement  in  the  objective  function  was  too  small.  As 
Algorithm  3.2  requires  a  finitie  upper  bound  (for  calculation  of  the 
step  size  in  step  2)  we  used  an  initial  value  of  UB  1.1*LB. 

Thereafter  for  UB  we  used  the  best  upper  bound  generated  so  far.  The 
parameters  and  tolerance  used  in  all  our  testing  were  these: 
c  =  .90 
NOLB  =  5 

NOUB  =  5 


Pricing  Tolerance  =  1.E-2 


▼ 
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4.2  Description  of  the  Test  Problems 

The  multicommodity  network  problem  generator,  MNETGN,  was  used  to 
create  22  random  test  problems.  We  modified  the  MNETGN  output  so  that 
every  arc  appeared  in  every  commodity's  subproblem  by  adding  arcs  with 
upper  bounds  of  zero  where  necessary.  The  test  problem  ranged  in  size 
from  22  to  754  nodes  and  from  53  to  1,102  arcs  with  from  0  to  599 
linking  constraints  and  from  3  to  20  commodities.  The  equivalent  LP 
sizes  are  between  232  and  8,904  rows  and  between  470  and  12,111  columns. 
The  number  of  linking  constraints  corresponds  to  a  wide  variety  of 
problems  from  pure  network  problems  (no  linking  constraints)  to  problems 
in  which  over  755®  of  the  arcs  are  included  in  linking  constraints. 

Problem  15  was  provided  by  Lt.  Col.  Dennis  McLain,  the  Assistant 
Director  of  Operations  Research  at  the  Military  Airlift  Command  located 
at  Scott  Air  Force  Base. 

4.3  Summary  of  Computational  Results 

All  the  testing  (except  for  Problems  15,  21,  and  23)  was  done  on  a 
CDC  6600  at  Southern  Methodist  University,  using  the  FTN  compiler  with 
the  optimization  feature  enabled.  Except  for  Problems  7  and  23,  a 
guaranteed  e -optimum  was  obtained  for  each  problem  with  t  >_  90S!. 

Problem  7  experienced  convergence  difficulties  when  run  using  EVAC. 
Problem  8  was  created  from  Problem  7  by  increasing  the  linking 
constraint  bounds  by  105®.  As  indicated  in  Table  4.1,  this  slight 
modification  enabled  EVAC  to  solve  the  problem  easily.  We  limited  the 
number  of  lower  bound  iterations  and  upper  bounds  iterations  to  100, 
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even  though  Problem  7  had  not  achieved  90%  optimality  by  that  point. 
Because  of  this  the  solution  times  for  Problem  7  are  given  in  Table  4.1 
but  are  not  included  in  the  summary  data. 

Problem  23  was  created  to  allow  us  to  test  EVAC  on  a  relatively 
large  problem.  This  problem  (with  8,904  LP  rows  and  12,111  LP  columns) 
was  too  large  for  MCNf  to  solve  in  the  available  memory,  so  we  were  not 
able  to  compare  solution  times  for  the  two  codes  on  this  problem.  In 
addition,  due  to  the  memory  limitations  on  the  CDC  6600,  we  were  forced 
to  use  a  CDC  203  to  test  Problem  23.  For  this  reason  the  times  for 
Problem  23  are  included  in  Tables  4.1  and  4.2,  but  are  not  included  in 
the  totals  and  summary  information.  Since  the  testing  on  the  CDC  203 
involved  a  real-dollar  expense,  we  were  satisfied  to  stop  when  a  lb% 
optimum  was  attained.  The  test  runs  for  Problems  13  and  21  were  made  on 
a  CDC  Cyber  73.  But  since  both  the  EVAC  and  MCNF  runs  for  these 
problems  were  made  on  the  Cyber  73,  the  totals  and  summary  data  include 
the  times  for  Problems  15  and  21. 

Details  of  the  test  problems  are  given  in  Table  4.1.  The  times  are 
in  CPU  seconds  and  exclude  the  time  required  to  input  the  problem  data 
and  print  the  solution  reports.  Table  4.1  also  presents  a  comparison  of 
the  times  required  for  MCNF  and  EVAC  to  solve  each  problem.  In  order  to 
present  a  meaningful  comparison  of  the  solution  times  for  MCNF  and  EVAC, 
we  also  present  the  solution  times  for  EVAC  exclusive  of  the  extra  1/0 
required  to  maintain  the  costs,  bounds,  and  old  bases  for  the  sub- 
problems  on  peripheral  storage.  Since  MCNF  maintains  all  this  informa¬ 
tion  in  main  memory,  this  seems  to  be  the  most  reasonable  way  of 
comparing  timing  statistics.  The  column  titled  "Guaranteed  %  Optimal" 
gives  the  best  lower  bound  generated  by  EVAC  as  a  percent  of  the  best 


upper  bound  generated  by  EVAC.  The  column  titled  "Actual  %  Optimal" 
presents  the  actual  optimal  objective  (as  obtained  by  MCNF)  as  a  percent 
of  the  best  upper  bound  generated  by  EVAC. 

Table  4.2  provides  the  details  of  the  times  required  by  EVAC  to 
perform  various  steps  of  the  algorithm.  The  column  titled  "%  of  Time  in 
Other"  for  the  lower  bound  computations  shows  the  time  required  for  such 
activities  as  computing  the  Lagrange  multipliers,  updating  the  unit 
costs  to  reflect  these  changes,  computing  the  resulting  dual  variables, 
and  various  bookkeeping  activities.  The  corresponding  column  for  upper 
bound  computations  reflects  such  activities  as  calculating  the  dual 
variables,  testing  the  termination  criteria,  and  various  other  short 
computations. 

Table  4.3  summarizes  the  time  comparisons  graphically.  The 
problems  are  grouped  by  number  of  commodities,  as  they  are  in  Tables  4.1 
and  4.2. 

4,4  Analysis  of  Results 

It  seems  clear  from  Tables  4.1  and  4.3  that  EVAC  severely 
dominates  MCNF  whenever  the  number  of  commodities  is  small.  This  is  due 
to  the  fact  that,  for  EVAC,  quite  a  bit  of  additional  overhead  is 
involved  in  alternating  between  commodities.  This  overhead  is  not  just 
a  result  of  I/O,  although  that  is  a  great  deal  of  it,  but  is  also  due  to 
the  set-up  time  required  for  activities  such  as  constructing  a  new 
feasible  basis  from  an  old  basis  and  calculating  the  resulting  dual 
variables.  MCNF,  on  the  other  hand,  is  primarily  driven  by  the  number 
of  binding  linking  constraints  in  the  optimal  solution.  This  is  because 
MCNF  seeks  an  exact  optimum. 
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Letting  T(EVAC)  denote  the  average  time  required  by  EVAC 
(exclusive  of  I/O),  and  T(MCNF)  denote  the  average  time  required  by 
MCNf,  we  can  express  the  following  relationships: 

For  the  3-commodity  test  problems, 

T(EVAC)  =  .354  *  T(MCNF). 

For  the  4-commodity  test  problems, 

T(EVAC)  =  .469  *  T(MCNF). 

For  the  5-commodity  test  problems, 

T(EVAC)  =  .666  *  T(MCNF). 

And  for  the  test  problems  with  6  or  more  commodities, 

T(EVAC)  =  .975  *  T(MCNF). 

It  should  also  be  noted  that  EVAC  is  capable  of  solving  larger 
problems  than  MCNF.  This  is  due  to  the  fact  that  EVAC  stores  only  one 
copy  of  the  network  defining  data  in  main  memory,  where  MCNF  requires 
one  copy  for  each  commodity.  Also,  EVAC  maintains  in  main  memory  the 
current  basis,  cost  and  bound  data  for  only  one  commodity  at  a  time. 
Thus,  for  a  K-commodity  problem,  EVAC  uses  on  the  order  of  1/K  the  main 
memory  required  by  MCNF. 

Note  that  the  entries  in  the  "Guaranteed  %  Optimal"  and  "Actual  % 
Optimal"  columns  of  Table  4.1  are  quite  close.  This  indicates  that  the 
sequence  of  lower  bounds  converged  to  values  very  near  optimality.  In 
addition,  from  Table  4.2,  we  see  that  the  lower  bound  iterations  are 
typically  less  time  consuming  than  the  upper  bound  iterations. 

It  is  worth  observing  that  EVAC  was  designed  for  very  large 
problems  which  would  never  be  solved  to  optimality.  Even  if  a  problem 
does  not  converge  to  within  the  requested  tolerance  in  a  prescribed 
number  of  iterations,  EVAC  always  provides  a  feasible  solution  which  is 
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a  guaranteed  e-optimum  for  somee>0.  In  contrast,  MCNF  provides  only  an 
upper  bound  on  the  optimum  objective  value,  with  no  indication  of  how 
close  it  is  to  optimality  until  an  exact  optimum  is  actually  attained. 

We  conclude  that  EVAC  works  extremely  well  in  obtaining  a 
guaranteed  e-optimum  for  the  multicommodity  network  flow  problem.  While 
it  is  not  as  "robust"  as  the  simplex-based  MCNF,  it  is  a  good  choice  for 
the  class  of  problems  for  which  it  was  developed,  the  very  large 
casualty  evacuation  models. 


TABLE  4.1 

DESCRIPTION  Of  THE  TEST  PROBLEMS  AND  SUMMARY  COMPARISON  Of  SOLUTION  TIMES  FOR  EVAC  AND  MCNF 


TABLE  A. 2 

DETAILED  TIMING  STATISTICS  FOR  EVAC  RUNS 
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TABLE  4.3 

GRAPHICAL  COMPARISON  OF  EVAC  AND  MCNF  SOLUTION  TIMES 
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CHAPTER  V 

SUMMARY  AND  CONCLUSIONS 


This  chapter  presents  a  summary  of  the  results  reported  in 
Chapter  IV  and  shares  conclusions  regarding  the  relative  effectiveness 
of  our  technique.  It  also  includes  ideas  for  further  investigation  in 
the  area. 

5.1  Summary  and  Conclusions 

Algorithm  3.6  describes  our  technique  for  finding  an  e-optimal 
solution  for  the  multicommodity  network  flow  problem.  Our  technique 
differs  from  other  approaches  to  the  problem  in  that,  rather  than 
solving  the  multicommodity  problem  directly,  we  compute  sequences  of 
lower  and  upper  bounds  on  the  optimal  objective  function  value, 
terminating  when  the  bounds  are  within  a  prescribed  tolerance.  Both 
the  lower  and  upper  bound  algorithms  use  a  subgradient  optimization 
technique  and  both  decompose  on  commodities  so  that  only  a  single 
commodity  minimum  cost  network  flow  optimizer  is  required.  At  each 
iteration  of  the  lower  bound  routine  (Algorithm  3.2),  an  initial  basis 
is  generated  from  the  previous  optimal  basis  by  modifying  the  costs  to 
correspond  to  the  new  Lagrange  multipliers,  and  updating  the  dual 
variables.  At  each  iteration  of  the  upper  bound  routine  (Algorithm 
3.5),  an  initial  basis  is  constructed  from  the  previous  optimal  basis 
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using  the  rules  described  in  [1]  to  restore  feasibility  (if 
necessary)  after  changing  the  bounds  to  correspond  to  the  new 
allocations. 

The  subgradients  for  the  lower  bounds  are  computed  to  be  the  sum 
of  the  flows  on  the  mutually  constrained  arcs  minus  the  associated 
mutual  arc  capacities.  For  the  upper  bounds,  subgradients  are 
computed  using  the  dual  variables  obtained  when  solving  the  single 
commodity  network  problems. 

□ur  computational  work  included  solving  each  one  of  23  problems 
twice;  once  using  MCNF,  a  primal  partitioning  code,  and  once  using 
EVAC,  our  implementation  of  Algorithm  3.6.  On  the  average  EVAC 
required  only  65SS  of  the  time  required  by  MCNF  (ignoring  J/0).  EVAC's 
performance  was  far  superior  on  the  problems  with  fewer  commodities 
and  was  not  as  impressive  on  the  problems  involving  many  commodities. 
In  addition  EVAC  required  on  the  order  1/K  the  amount  of  main  memory 
as  MCNF  for  a  K-commodity  problem. 

5.2  Areas  for  Future  Investigation 

Algorithm  3.6  involves  two  more  or  less  independent  processes. 
That  is,  there  is  no  reason  why  the  lower  bound  generator  (Algorithm 
3.2)  and  the  upper  bound  generator  (Algorithm  3.5)  could  not  proceed 
independently,  stopping  now  and  then  to  exchange  their  best  bounds  and 
test  for  optimality.  Hence  it  appears  that  this  procedure  is 
well-suited  to  exploit  the  benefits  of  a  parallel  processing 
environment.  In  addition  to  the  partitioning  of  the  technique  into 
two  separate  procedures,  within  each  of  these  procedures  the 
decomposition  by  commodities  could  take  advantage  of  a  parallel 
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processing  scheme  as  well.  It  would  seem  reasonable  to  expect  such  a 
scheme  to  speed  up  the  execution  time  considerably,  especially  when 
solving  a  very  large  problem. 

There  is  also  room  for  additional  experimentation  with  the  step 
sizes,  specifically  with  the  multipliers  on  the  step  sizes.  Perhaps  a 
scheme  in  which  the  multipliers  were  allowed  to  be  reset  to  their 
starting  values  a  finite  number  of  times  would  speed  up  convergence. 
One  might  reset  these  multipliers  whenever  the  improvement  in  the 
sequence  of  upper  (lower)  bounds  fell  below  some  tolerance.  This 
would  have  the  effect  of  restarting  the  algorithm  at  that  point,  but 
with  a  far  better  "starting  solution". 

In  addition  this  problem  has  a  multiperiod  structure.  Since  the 
network  is  replicated  for  60  one  day  time  periods,  it  might  be 
advantageous  to  exploit  this  structure  using  a  forward  simplex 
approach. 
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ABSTRACT 


This  paper  presents  a  new  algorithm  to  solve  a  network  problem 
with  equal  flow  side  constraints.  The  proposed  solution  technique  is 
motivated  by  the  desire  to  exploit  the  special  structure  of  the  side 
constraints  and  to  maintain  as  much  of  the  characteristics  of  pure 
network  problems  as  possible.  Not  only  has  specialized  software  for  the 
efficient  solution  of  pure  networks  been  developed,  but  the  same  compu¬ 
tational  efficacies  lend  themselves  to  the  solution  of  sequences  of 
minimum  cost  network  flow  problems  by  using  reoptimization  procedures. 

Our  solution  technique  for  the  equal  flow  problem  consists  of  solving 
two  sequences  of  pure  network  problems.  One  sequence  yields  tighter 
lower  bounds  on  the  optimal  value  by  considering  the  Lagrangean  relax¬ 
ation  of  the  equal  flow  problem  in  which  the  side  constraints  are  dualized. 
The  second  sequence  yields  upper  bounds  on  the  optimal  value  for  the 
problem  and  maintains  a  feasible  solution  at  all  times.  This  sequence  is 
obtained  by  considering  a  reformulation  of  the  equal  flow  problem  based 
on  parametric  changes  in  the  requirements  vector.  The  procedure  has  the 
added  attractive  feature  that  it  provides  a  feasible  solution  which  is 
known  to  be  within  a  percentage  of  the  optimal  at  all  times.  As  6uch, 
the  algorithm  terminates  when  a  solution  with  a  prespecified  tolerance 
on  the  objective  function  value  is  obtained.  On  NETGEN  problems, 
using  the  first  150  arcs  to  form  75  equal  flow  side  constraints,  we  found 
that  the  new  algorithm  is  approximately  3  times  faster  than  existing 
techniques  and  requires  only  50%  of  the  storage. 
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I.  INTRODUCTION 


This  paper  presents  a  new  technique  to  solve  the  equal  flow  problem. 
This  problem  is  easily  conceptualized  as  a  minimal  cost  network  flow 
problem  with  additional  constraints  on  certain  pairs  of  arcs.  Specifically, 
given  pairs  of  arcs  are  required  to  take  on  the  same  value.  Applications 
of  this  model  include  crew  scheduling  [6],  estimating  driver  costs  for 
transit  operations  [28],  and  the  two  duty  period  scheduling  problem  [25]. 

The  equal  flow  problem  may  be  solved  using  a  specialization  of  the  simplex 
method  for  networks  with  side  constraints.  However,  by  exploiting  the 
special  structure  of  the  side  constraints,  we  have  developed  a  new  algo¬ 
rithm  which  results  in  a  decrease  in  both  computer  storage  and  compu¬ 
tation  time. 

It  is  well  documented  that  pure  network  problems  can  be  solved  from 
fifty  to  one  hundred  tiroes  faster  using  specialized  primal  simplex  soft¬ 
ware  as  compared  to  general  linear  programming  systems.  Motivated  by  this 
great  advantage,  our  procedure  solves  the  equal  flow  problem  as  a  sequence 
of  pure  network  problems  and  totally  eliminates  the  need  to  deal  with  a 
basis  matrix. 

1.1  Problem  Description 

The  equal  flow  problem  is  defined  on  a  network  represented  by  an  (m,n) 
node-arc  Incidence  matrix,  A,  in  which  K  pairs  of  arcs  are  identified  and 
required  to  have  equal  flow.  Mathematically,  this  is  expressed  as: 


Minimize  cx 


s.t.  Ax  *  b 

Xk  “  *k+K  *  k  -  1,  ....  K 

0  <  x  <  u 

where,  c  Is  a  1  x  n  vector  of  unit  costs,  b  is  an  m  x  1  vector  of  node 
requirements,  0  is  an  n  x  1  vector  of  zeroes,  x  is  an  n  x  1  vector  of 
decision  variables,  and  u  is  an  n  x  1  vector  of  upper  bounds.  The  above 
definition,  henceforth  referred  to  as  PI,  assumes  that  the  first  2K  arcs 
appear  in  the  equal  flow  constraints.  This  assumption  is  in  no  way 
restrictive  since,  by  rearranging  the  order  of  the  arcs,  any  equal  flow 
problem  with  K  pairs  can  be  expressed  in  the  above  form.  Note  that  the  K 
pairs  of  arcs  are  mutually  exclusive,  i.e.  an  arc  appears  in  at  most 
one  side  constraint. 

1.2  Survey  of  Related  Literature 

In  1961,  Chames  and  Cooper  (7]  presented  a  specialized  algorithm 
for  the  model: 


Minimize 

cx 

Set* 

Ax  *  b 

Cx  -  d 

x  >  0, 

where  A  and  C  are  some  general  matrices  but  A  has  some  favored  structure. 
Their  algorithm,  called  the  double  reverse  method,  takes  advantage  of 
the  special  structure  of  the  matrix  A.  Variations  of  this  algorithm  may 
be  found  in  [2,  8,  10,  15,  19,  24].  Specializations  for  multicommodity 
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problems  may  be  found  in  [14,  20,  21]. 


In  1980,  Shepardson  and  Marsten  [25]  showed  that  the  two  duty  period 
scheduling  problem  can  be  reformulated  as  a  single  duty  period  scheduling 
problem  with  e^ual  flow  side  constraints.  They  obtain  a  Lagrangean  dual 
for  this  equal  flow  problem,  by  dualizing  with  respect  to  the  equal  flow 
side  constraints.  This  Lagrangean  dual  is  maximized  using  the  subgradient 
optimization  technique.  In  1984,  Turnquist  and  Malandraki  [28]  modeled  the 
problem  of  estimating  driver  costs  for  transit  operations  as  an  integer 
equal  flow  problem.  They  obtain  a  Lagrangean  dual  for  their  problem,  by 
dualizing  with  respect  to  the  side  constraints.  Their  algorithm  is  a 
slight  modification  of  the  subgradient  optimization  technique.  They 
perform  a  line  search  between  two  successive  solutions  obtained  during 
the  subgradient  optimization  process. 

Beck,  Lasdon,  and  Engquist  [5]  transformed  the  equal  flow  problem 
into  a  quadratic  programming  problem  which  has  a  penalty  for  violating 
the  equal  flow  constraints.  They  solved  this  nonlinear  programming 
problem  using  the  Fletcher-Reeves  conjugate  gradient  method  [9],  a 
successive  linear  programming  code  [13],  and  a  convex  simplex  code. 

If  the  penalty  is  sufficiently  large,  this  approach  is  guaranteed  to 
converge  to  the  optimal  solution  of  the  equal  flow  problem. 

1,3  Objective  of  the  Investigation 

The  objective  of  this  investigation  is  to  develop  and  computationally 
test  a  new  algorithm  for  the  equal  flow  problem.  This  algorithm  utilizes 
the  subgradient  optimization  technique  and  is  based  on  the  relaxation/ 
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restriction  procedure  proposed  by  Glover,  Glover,  and  Martinson  [11]  for 
a  generalized  network  model  with  special  side  constraints.  We  establish 
that  the  equal  flow  problem  may  be  solved  as  two  sequences  of  pure 
network  problems,  one  sequence  corresponds  to  computing  a  lower  bound 
while  the  other  corresponds  to  computing  an  upper  bound.  In  the  limit, 
both  bounds  will  converge  to  the  optimal  objective  value.  Our  implementation 
terminates  when  the  difference  between  the  bounds  is  within  a  prespecified 
tolerance. 

The  subgradient  optimization  technique  requires  the  computation  of 
6ubgradients,  choice  of  appropriate  step  sizes,  and  the  application  of 
a  projection  operation.  We  show  that  the  subgradients  for  the  upper 
bound  can  be  computed  using  the  optimal  dual  variables  obtained  by 
solving  pure  network  problems.  We  also  develop  theoretical  results 
that  yield  an  easy  implementation  of  the  projection  operation.  The 
step  sizes  selected  are  a  modification  of  the  ones  proposed  by  Polyak 
[23].  For  this  choice  of  step  sizes,  we  prove  that  our  algorithm  must 
necessarily  obtain  an  iterate  at  which  the  objective  value  is  arbitrari¬ 
ly  close  to  the  optimal  objective  value.  In  a  computational  study, 
comparing  our  code  with  a  code  that  is  designed  to  solve  network  problems 
with  side  constraints,  we  found  that  the  new  code  runs  approximately  3 
times  faster  and  requires  30%  less  core  storage. 
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II.  THE  SUBGRADIENT  ALGORITHM 


The  Subgradient  Algorithm  was  first  introduced  by  Shor  [27]  and  is 
a  general  procedure  for  solving  nonlinear  programming  problems.  It  may 
be  viewed  as  a  generalization  of  the  steepest  descent  (ascent)  method 
for  convex  (concave)  problems  in  which  the  gradient  may  not  exist 
everywhere.  The  subgradient  is  simply  substituted  in  place  of  the 
gradient  for  those  points  for  which  the  gradient  does  not  exist.  When 
this  occurs,  the  algorithm  may  move  to  a  point  with  objective  value 
worse  than  the  current  point.  Hence,  the  objective  function  does  not 
necessarily  improve  at  each  iteration  and  consequently  the  convergence 
results  of  Zangwill  [29]  do  not  apply.  Remarkably  though,  under  fairly 
minor  conditions  on  the  step  size,  convergence  can  be  guaranteed. 

Let  the  nonlinear  program  PO  be  given  by: 

Minimize  f(y) 
s.t.  y  c  G 

where  f  is  a  real  valued  function  that  is  convex  over  the  compact,  con¬ 
vex,  and  nonempty  set  G.  A  vector  n  will  be  called  a  subgradient  of  f 
at  y  if  f(y)  -  f(y)  >  r)(y  -  y)  for  all  y  E  G,  For  any  y  £  G,  we  denote 
the  set  of  all  subgradients  of  f  at  y  by  3f(y).  The  subgradient  algorithm 
makes  use  of  an  operation  called  the  projection  operation.  The  projection 
of  a  point  x  onto  G,  denoted  by  P[x],  is  defined  to  be  the  unique  point  y  e  G 
that  is  nearest  to  x  with  respect  to  the  Euclidean  norm.  Using  the  projec¬ 
tion  operation,  we  now  present  the  subgradient  algorithm  in  its  most 
general  form. 
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ALG  1  SUBGRAVJENT  OPTIMIZATION  ALGORITHM 


Step  0  (Initialization) 

Let  y^  be  any  element  of  G,  select  a  set  of  step  sizes 

s0,sl,s2****'  an<*  set  *  ■*" 

Step  1  (Find  Subgradient) 

Let  e  3f(y^).  If  *  0,  then  terminate  with  y  optimal. 

Step  2  (Move  to  New  Point) 

Set  y±+1  *■  Pfy^  -  s^tl],  set  i  ■*-  i  +  1  and  return  to  step  1. 

Various  proposals  have  been  offered  for  the  selection  of  the  step 
sizes.  Three  general  schema  which  have  been  suggested  are: 


i) 


V 


ii) 


A  (f(y  >-f*) 

iii)  s,  «  — - ^ -  ,  0  <  A.  <  2, 

1  IKII2  1 


where  f*  is  the  optimal  value  of  f  over  G.  If  the  constants,  A^s, 
satisfy  the  following  conditions: 

1  *11  **  11®  1*  *  0;  and  £  A.  *  •*>, 

i-*  *  1 

then  the  convergence  of  the  algorithm  is  guaranteed  using  <i)  or  (ii) 
(see  Coffin  [12],  Helgason  [17],  Kennington  and  Helgason  [21]).  For  the 
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upper  bounds,  we  use  a  modification  of  the  third  step  size.  The 
following  result  is  available  for  this  scheme. 

Proposition  1  (Polyak  [23]) 

Let  f  be  a  real  valued  convex  function  over  the  compact,  convex, 
and  nonempty  set  G.  Also,  let  f*  be  the  minimum  of  f  and  ||n^||  <  C  for 
all  i  and  some  constant  C.  Then  there  exists  a  y*  e  G  with  f(y*)  •+  f*, 
if  scheme  (iii)  is  used. 

Note  that  in  (iii),  f*  is  the  optimal  value  of  f  over  G.  Since  the 
optimal  objective  is  unknown  before  solving  the  problem,  we  use  a  lower 
bound  on  f*  in  our  implementation. 
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III.  THE  LOWER  BOUND 


Recall  that  the  equal  flow  problem,  which  we  denote  by  PI,  is 
given  by: 

Minimize  cx 
s.t.  Ax  -  b 

*k  *  *K+k  ’  k  “  1*  — »K 
0  <  x  <  u. 

In  our  algorithm  for  PI,  lower  bounds  on  the  optimal  objective  of  PI  are 
used  for  step  sizes  and  for  termination.  In  this  section,  we  describe  a 
procedure  to  obtain  these  lower  bounds . 

Consider  the  following  Lagrangean  dual  for  PI,  which  we  shall 

refer  to  as  Dl: 

Maximize  h(v)  ,  where  w  *  [w^ . wR]  C  R v,  and 

K 

h(w)  -  Min{cx  +  l  “  XK+^ :  Ax  E  b,  0  <  x  <  u}. 

k=l 

Proposition  2  (Shetty  [26]) 

Let  x  be  a  feasible  solution  to  PI  and  let  w  *  [w^,...,w^]  be 
a  feasible  solution  to  Dl.  Then  cx  >  h(w). 

Proposition  3  (Bazaraa  and  Shetty  [4]) 

If  PI  has  a  minimum,  then  the  optimal  objectives  for  PI  and  Dl 

are  equal. 

As  a  consequence  of  Propositions  2  and  3,  we  may  solve  Dl  to 
obtain  a  lower  bound.  We  will  now  show  that  Dl  may  be  solved  using  the 


subgradient  optimization  technique  for  concave  functions.  This  tech¬ 
nique  is  similar  to  ALG  1  with  a  modification.  Let  Pq»P^,P2»*** 
denote  a  sequence  of  step  sizes  and  let  di  e  Sh^).  Then  step  2  is 
replaced  by: 

Step  2  (Move  to  New  Point) 

Set  +  P^d^,  set  *  and  return  to  step  1. 

To  use  this  algorithm  h(w)  oust  be  concave,  and  we  need  a  means 
of  generating  subgradients.  These  two  results  follow: 

Proposition  4  (Shetty  [261) 

ir 

The  real  valued  function  h  is  concave  over  R  . 

Proposition  5  (Shetty  [26]) 

For  a  given  w,  let  x  be  an  optimal  solution  to 
K 

Min{cx  +  J  w 
k=l 

Then  d  *=  [(x^  -  XK+1>  » •  •  • «  (*K  “  X2K^  iS  3  subRradlent  of  h  at  *’• 

We  used  scheme  (i)  for  step  sizes.  Let  UBND  denote  an  upper  bound  and 
assume  that  the  optimal  objective  value  is  positive.  Our  algorithm  for 
obtaining  lower  bounds  is  presented  below: 

A LG  2  LOWER  BOUND  ALGORITHM 

Step  1  (initialization) 

Initialize  UBND,  step  size  p,  and  tolerance  e. 

Set  w  *■  0. 


k(xk  "  *K+k):  ^  =  b»  0  i  x  t  u)« 
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Step  2  (Find  Subgradient) 

Let  x  *  [x  ,...,xnl  solve 

h(v)  «  Min{cx  +  f  vk(xk  -  x^) :  Ax  «  b,  0  <  x  <  u}. 
k=l 

Set  LEND  h(w). 

If  (UBND  -  LBND)  <  e(UBND)  then  terminate; 
otherwise,  set  d  *•  10^  -  x^) » •  •  • » (*£  ~  *2K^' 

Step  3  (Move  to  a  New  Point) 

3a.  Set  w  w  +  pd,  set  p  p/2. 

3b#  Go  to  2# 
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IV.  THE  UPPER  BOUND 


An  alternate  formulation  of  PI,  which  will  be  referred  to  as  P2, 
is  as  follows: 


Minimize  g(y) 
s.t.  y  c  S 

where  for  any  vector  y  ■  ly^,..,,y  ], 

g(y)  ■  Min{cx:  Ax  =  b,  0  <  x  <  u,  ■  xK+k  *  ^or  ^ 

and 


S  «  {y:  0  <  y^  <  min(uk,uK+k>  for  all  k}. 

Clearly,  PI  and  P2  are  equivalent.  That  is,  given  an  optimum  for  one, 
we  can  construct  an  optimum  for  the  other.  We  will  now  show  that  P2 
is  a  special  case  of  the  nonlinear  program  PO  and  may  be  solved  using 
the  Subgradient  Optimization  Algorithm,  ALG  1. 

Proposition  6  (Shetty  [26]) 

The  real  valued  function  g  is  piece-wise  linear  convex  over  the 
compact,  convex  and  nonempty  set  S. 

To  apply  the  subgradient  algorithm,  we  need  a  procedure  for  obtaining 
a  subgradient  of  g  at  a  point  y.  The  following  proposition  shows  that 
the  dual  variables  may  be  used  to  construct  a  subgradient. 


Proposition  7  (Shetty  126]) 

Let  <*»VVK+1,,**,VV2K’ 


y)  be  the  optimal  dual  variables  for 
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Minimize  cx 


*K+1  "  yl 


(VK+1J 


2K  y2K 


(v2K> 


0  <  x  <  u  (p). 

Then  n  *  (v.^  +  vK+j»**«»vK  +  V2K)  is  a  subgradient  of  g  at  y  “  (y1,...,yJC). 

As  a  result  of  Proposition  7,  a  subgradient  at  any  given  point 
y  *  (y^,...,yK>  e  s  required  in  our  specialization  of  ALG  1  can  be 
obtained  by  solving  Min{cx:  Ax  *  b,  x^  *  ^l’***,x2K  *  yK’  ®  -  *  -  u^’ 
we  shall  refer  to  as  P3.  After  substituting  x^  *  y^,...,X2K  *  yR*  in 
Ax  *  b,  we  obtain  a  pure  network  problem,  which  we  shall  refer  to  as 
P4  and  is  given  below: 

AA 

Minimize  cx 

AA  /\ 

s.t.  Ax  -  b 


0  <  x  <  u. 

To  apply  ALG  1,  we  need  a  procedure  for  constructing  the  optimal  dual 
variables  (vi»V£+i» • • • ,vk*v2K^  ^or  ^  *rom  OP1*®*!  dual  variables 
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for  P4.  Suppose  the  arc  corresponding  to  j  has  "From"  node  j  and  "To" 
node  J2*  T^at  is»  arc  J  *s  the  ordered  pair  (J^ » ^2^ ”  Then  we  define 
FROM(j)  *  and  TO(j)  «  Using  this  notation,  the  following  propo¬ 

sition  gives  the  required  formulae: 

Proposition  8  (Shetty  [26]) 

A 

Let  tt  be  the  vector  of  optimal  dual  variables  for  P4.  Then 

[IT,  vi*vk+i»**‘»vk»v2k]  With  vj  "  “^FROM(j)  +  *TO(j)+  C  j  ’ 
are  optimal  duals  for  P3. 

We  now  present  two  propositions  that  justify  the  projection  routine 
used  for  the  upper  bound.  The  proofs  may  be  found  in  Kennington  and 
Helgason  [21]  and  Shetty  [26], 

Proposition  9 

Let  S  be  a  nonempty,  convex  set  and  y  t  S.  Then  y*  e  S  is  a 
projection  of  y  on  to  S  if  ( y  -  y*) (y  -  y*)  <  0  for  all  y  e  S. 

Proposition  10 

Let  y  -  (ylty2f»,yK)  e  RK  with 

yk  <  0  for  k  ■  1,...,L 

0  :  yt  5  “k.  \  ■  ®ln(uk,  uK+k)  for  k  -  I-+1 . L+M 


yk  >  uk 


for  k  -  L+M+1,...,K 
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where  L,M  are  Integers  and  0  <  L,M  <  K.  Then 


P(y)  -  y*  -  (0, . . . ,o,?L+i, . . .  ,yL+M.uL4M+1 . V  18 

A 

a  projection  of  y  on  S* 

Following  a  description  of  the  terminology  used,  our  algorithm  for 
obtaining  an  upper  bound  for  PI  is  presented  below.  Let  RFREQ  denote  the 
frequency  at  which  the  constant  X  in  step  size  (iii)  is  reset  to  its  initial 
value,  Xq,  LBND  denote  a  lover  bound  on  PI,  UBND  denote  an  upper  bound  on 
PI,  P  denote  the  projection  routine  described  in  Proposition  10,  e  denote 
the  termination  tolerance  and  Q  denote  the  iteration  count  for  the 
upper  bound . 


A LG  3  UPPER  B0UMP  ALGORITHM 

Step  1  (Initialization) 

Choose  yeS. 

Initialize  LBND,  RFREQ,  e,  and  X^. 

Set  u  (min(u1,uK+1>  , . . .  ,min(uK,u 

Set  Q  0,  set  ^  «-  XQ  for  k»l,...,K. 

Step  2  (Find  Subgradient  and  Step  Size) 

For  allocation  y,  let  x  and  it,  respectively,  be  the  vectors  of 

optimal  primal  and  dual  variables  for 

Min{cx:  ££■=£,  0  <  x  <  u}.  Construct  x  from  x,y. 


Set  UBND  +■  cx 


If  (UBND  -  LBND)  <  e(UBND),  then  terminate  with  x  optimal; 
otherwise, 

set  Vj  +  ^(j)  +  C j  *  j  1 . 2K’ 

set  n  -  (vx  +  vK+1,...,vKfv2K). 

If  Q  *  RFREQ,  then  set  Q  0,  6et  X  Xq,  set  «*-  Xq 

for  k  *  and  go  to  3; 

otherwise , 

compute  X^  such  that  0  <  y^  -  X^r^  5  ujc*  k  * 
set  min{I^/2,Xk},  k  ■ 

set  X  min{X^,k*l,2, . . .  ,K}. 

Step  3  (Move  to  New  Point) 

3e.  Set  y  *■  Ply  -  A  (l)BK°  LBf))  n] ,  set  Q  -e  Q+l. 

II  tl  ||  2 

3b.  Go  to  2. 

Note  that  the  step  size  (iii)  presented  before  may  be  rewritten 
for  our  function  g  as  follows: 

^i(g(yi)  -  8*) 

s,  -  -  ,  o  <  X  <  2, 

II  nA  II2 

where  r|±  C  SgCy^  and  g*  is  the  optimal  value  of  g. 
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In  our  implementation,  we  use  g,  a  lower  bound  on  g,  in  place  of  g*. 

The  following  propositions  demonstrate  that  for  g’s  close  to  g*,  our 
procedure  must  necessarily  obtain  an  iterate  at  which  g  is  arbitrarily 
close  to  g*. 

Proposition  11  (Kennington  and  Helgason  [21]) 

if  ni  +  0, 

II  ~  yj II 2 1  II  y±  ~  y II 2  +  si II  II 2  +  2sir>i (y  ■  yA)  f °r  any  yes 

and  step  size  s^. 

Proposition  12 

Let  g*  be  the  optimal  value  of  g,  and  also  let 

i)  a  g*  <  g  <  g*,  0  <  a  <  1, 

ii)  s1  -  X1(g(y1)  -  I)/  ||  ndH  2  »  and 

iii)  0  <  e  <  ^<  B  <  2  for  all  i. 

If  there  is  a  constant  C  such  that  j|  Tl^[|  <  C  for  all  i,  then  there 

2  B 

exists  some  i  such  that  g(yA)  <  M6  +  g*  (■  - jqp  a)  for  any  6  >  0 

and  for  some  constant  M. 

Proof 

First,  we  assert  that  there  is  some  i  such  that 

A l  i  2X±  g*  d 

g(y  )  <  - - - - - - -  ,  for  any  6  >  0. 

1  (X±  -  2Ai)  (A‘  -  2A1)  (A‘  -  2Aa) 
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Suppose  that  for  all  i. 


x2  S 


2Xi  g* 


g<yi)  >  —2 


,  vhere 


a[  -  2Xi)  (X‘  -  2X±)  <Xi  -  2\±) 


6  >  0  is  given.  Let  y*  e  S  be  an  optimal  point.  By  Proposition  11, 

2X1(g(y1)-I) 


x^<g(yi)-i)2 


yi41-y*l!  5  «  V»*“  + 


II  nt 


Since  t\.  c  9g(yi),  -  y±)  5  8*  ' 


Thus, 


X^gty^-g)2  2X1(g(yi)-g)(g*  -  g(yi)) 


yi+l‘y 


-v*  r  < 


y±  -  y*  II  + 


li  ni 


II  ni 


yi  ~  y*  II2  + 


g) 


g(yi)(A±  "  2V~Xi  *  +  2Xi  8* 


ii  »»«  ir 


<g(yj  ~  g)<5 


y±  ~  y*  l!2  - 


II  ^  II ' 


-  y*  II2  - 


(g*  -  i)<5 


(i) 
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Me  can  choose  an  integer  N  large  enough  that 


c2  II  y*  -  y*  I!2 


(g*  “1)6 


<  N. 


Adding  together  the  inequalities  obtained  from  (1)  by  letting  i  take 
on  all  values  from  1  to  N,  we  obtain 

N(g*  -  1)6 


yN  ■  y*  II2  <  II  y±  -  y*  II2 


a  contradiction.  This  justifies  our  assertion. 
By  simplifying  our  assertion  further  we  get, 

6 


<  0, 


g(y1)  < 


8  K 


2g* 


(2A  -  A  ;  (2  -  A.)  (2  -  A.) 

ii  i  i 


a  g**.,  2g* 

- i -  +  - 


(2  X±  -  Ap  (2  -  At)  (2  -  A1) 


(2Xi  '  V 


<2Ai  -  A J) 


+  g* 


2  -  a  A, 


2  -  K 


+  g*(l  + 


(2  - 


(1  -  a)) 


6  0 
< - -  +  g*(1  + -  (1  _  a)) 


(2At  -  A‘) 


2-B 
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&  2  3 

■  - s —  +  8*( - °) 

(2Ai  -  A*)  2-3  2-3 

2  3 

<  M6  +  g*( - a),  where  M  is  a  constant  less  than 

j_  "  2-3  2-3 

2  * 

(2e-0 

This  completes  the  proof  of  Proposition  12. 
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V.  THE  ALGORITHM 


In  this  section,  we  present  our  new  algorithm  for  solving  the  equal 
flow  problem.  Let  ITERL  denote  the  number  of  iterations  spent  in  step  la 
in  computing  the  lower  bound  before  returning  to  the  upper  bcund,  and 
ITERU  denote  the  number  of  Iterations  spent  in  computing  the  upper  bound 
before  returning  to  the  lower  bound.  Also,  let  T  denote  the  iteration 
count  for  the  lower  bound,  R  denote  the  iteration  count  for  upper  bound, 
and  Pq  denote  the  initial  step  size  for  the  lower  bound. 

ALG  4  SUBGMMEUT  OPTIMIZATION  ALGORITHM 
TOR  THE  EQU PLOW  PROBLEM 

Step  0  (Initialization) 

Initialize  ITERL,  ITERU,  REREQ,  AQ,  pQ,  and  tolerance  e. 

Set  T  0,  set  Q  0,  set  R  *■  0,  set  w  *■  0,  and  set  IFLAG  *•  0. 

Set  UBND  *-  +  «*>,  set  LBND  *■  -  <*>,  set  p  ■*-  p^,  and  set  A^  *■  A^  for  k=l,...,K. 
Set  u  ■*-  (min(u1,uK+1),...,min(uK,u2K)). 

Step  1  (Compute  Lower  Bounds) 

la.  Call  ALG  2  (steps  2  and  3a). 

lb.  Set  T  <*-  T+l 

If  T  ITERL,  then  go  to  1. 

lc.  (Initialize  y) 

If  IFLAG  +  0,  then  go  to  2;  otherwise, 

set  y  (min(u1,(x1+xK+1)/2),...,min(uK,(xK+^2K)/2)). 
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Step  2  (Compute  Upper  Bounds) 


2a.  Set  T  «-  0,  aet  IFLAG  «-  1. 

2b.  Call  Alg  3  (steps  2  and  3a) 

2c.  Set  R  ■+-  R+l. 

If  R  <  1TERU,  then  go  to  2b. 

Step  3 

Set  R  «-  0. 

Set  p  •*-  pQ. 

Go  to  1. 

In  the  above  algorithm,  IFLAG  is  used  in  obtaining  a  starting  y  from  the 
solutions  in  la.  The  bases  used  in  steps  1  and  2  are  generated  from  the 
optimal  bases  obtained  in  the  previous  iterations. 
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VI.  COMPUTATIONAL  EXPERIMENTATION 

This  section  describes  the  computer  implementation,  EQFLO,  and 
testing  of  our  algorithm  for  the  equal  flov  problem.  The  algorithm 
was  tested  on  a  set  of  35  test  problems  randomly  generated  using 
NETGEN  [22].  Computation  times  are  compared  with  those  of  NETSIDE  [2], 
a  general  purpose  code  for  network  problems  with  side  constraints. 

Both  NETSIDE  and  EQFLO  are  written  in  standard  FORTRAN  for  an  incore 
implementation  and  have  not  been  tailored  to  either  the  machine  or 
FORTRAN  compiler  used  for  testing. 

6.1  Description  of  the  Computer  Codes 

NETSIDE  was  developed  by  Barr,  Farhangian,  and  Kennington  at 
Southern  Methodist  University,  Dallas,  Texas.  Designed  to  solve  network 
problems  with  side  constraints,  it  used  a  specialization  of  the  revised 
simplex  method  known  as  the  primal  partitioning  algorithm  [15],  The 
basis  inverse  is  maintained  as  a  rooted  spanning  tree  and  a  working  basis 
inverse  in  product  form.  The  reinversion  routine  is  a  modification  of  the 
work  of  Hellerman  and  Rarick  [18]  and  uses  the  "spike  swapping  theory" 
of  Helgason  and  Kennington  [16],  The  initial  working  basis  consists  of  a 
combination  of  artificial  and  slack  variables.  The  working  basis  is 
reinverted  every  50  iteiations.  The  pricing  routine  uses  a  candidate  list 
of  size  10  with  a  block  size  of  400.  Both  pricing  and  pivot  tolerance  are 
l.E-6. 

EQFLO  is  our  implement a i ton  of  ALG  4,  and  makes  use  of  MODFLO  [1]  to 
solve  pure  network  subproblems.  MODFLO  is  a  set  of  subroutines  which  may 
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be  used  to  solve  a  network  problem  as  well  as  reoptimize  after  problem 
data  changes.  Based  on  NETFLO  [21],  this  code  allows  the  user  to  change 
costs,  bounds  and/or  requirements  for  a  network  problem  and  reoptimize. 

The  tuning  parameters  used  in  all  runs  were  as  follows:  1TERL  ■  5, 

ITERU  -  10,  REFREQ  -  5,  XQ  -  0.75,  pQ  -  0.01,  and  e  -  0.1.  M0DFL0  [1] 
is  used  to  reoptimize  after  each  change  to  either  the  costs  or  right- 
hand-sides. 

6.2  The  Test  Problems 

The  program  NETGEN,  a  generator  for  large-scale  network  test 
problems,  was  used  to  generate  35  test  problems.  The  parameters  used 
to  generate  these  problems  are  described  in  Klingman,  Napier,  and  Stutz  [22]. 
The  test  problems  have  between  200  and  1500  nodes,  and  1500  and  7000  arcs. 

For  each  problem,  the  first  150  arcs  were  paired  to  form  equal  flow  sides 
constraints.  The  characteristics  of  these  test  problems  are  listed  in 
Table  1. 

Our  algorithm  requires  upper  bounds  on  all  equal  flow  arcs.  Though 
NETGEN  generates  bounds  on  some  of  these  arcs,  there  were  others  with  no 
upper  bounds.  We  set  the  maximum  of  all  supplies  and  demands  to  be  the 
upper  bounds  on  such  arcs.  These  bounds  were  acceptable  since  the  optimal 
solutions  obtained  for  our  pure  network  problems  were  the  same  as  the  ones 
listed  in  NETGEN.  Furthermore,  for  all  35  test  problems  the  first  150  arcs 
were  used  to  form  75  pairs  of  equal  flow  side  constraints.  We  were  unable 
to  experiment  with  more  than  75  pairs  due  to  a  core  storage  limitation  of 
301K.  NETS2DE  required  approximately  300K  octal  words  of  storage  for 
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problems  28  through  35  with  75  pairs  and  any  further  increase  in  the 
number  of  pairs  would  exceed  the  storage  limitation. 


Table  1  About  Here 


6.3  Computational  Results 

All  35  test  problems  were  solved  on  the  CDC  6600  at  Southern 
Methodist  University,  using  the  FTN  compiler  with  OPT  “2.  All  35 
problems  were  solved  twice  using  EQFLO;  once  with  the  same  step  size 
for  every  pair  of  equal  flow  arcs  and  the  second  time  with  different 
step  sizes  for  different  pairs.  While  using  EQFLO  to  solve  these 
problems,  ALG  4  was  followed  exactly  the  first  time,  whereas,  the 
computation  of  the  step  size  for  the  upper  bound  was  altered  the  second 
time.  The  modification  was  as  follows: 


Step  3  (Move  to  New  Point) 


3a.  Set  yR  •*-  P 


yk"  \ 


(UBND-LBND) 

a  \ 


11  n  ir 


,  k  m  1, • . • ,K, 


set  Q  Q+l. 


Mote  that  this  modification  results  in  different  step  sizes  for  different 
equal  flow  pairs.  The  details  of  all  runs  are  given  in  Tables  2  and  3. 
The  times  are  in  CPU  seconds  and  exclude  input  and  output. 

The  value  0.01  used  for  p^  worked  well  for  all  test  problems  except 
problem  number  11.  This  problem  experienced  difficulties  in  converging 
within  10%  of  the  optimal.  However,  the  problem  did  converge  within  10% 
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of  the  optimum  when  we  changed  Pq  to  values  between  3  and  10. 

Tables  2  and  3  About  Here 


The  computational  results  presented  in  Tables  2  and  3  are  summarized 
in  Table  4.  Letting  T(ALG)  denote  the  CPU  time  required  to  solve  the  35 
test  problems  using  code  ALG,  the  relationship  is  given  below: 

T(NETSIDE)  *  2.54  (EQFLO),  same  step  size, 

T(NETSIDE)  -  3.00  (EQFLO),  different  step  sizes. 

Note  that  EQFLO  performs  better  as  the  problem  size  increases.  Although 
EQFLO  was  slightly  slower  than  NETSIDE  on  problems  1  to  10,  its  performance 
increased  substantially  on  problems  11  to  35.  In  particular,  EQFLO  ran 
approximately  5  to  6  times  faster  than  NETSIDE  on  problems  28  through  35 
and  these  problems  are  fairly  large.  We  expect  EQFLO  to  perform  even 
better  on  much  larger  problems.  This  is  attributable  to  the  fact  that 
the  time  for  pricing  and  updating  increases  dramatically  for  NETSIDE 
with  an  increase  in  the  size  of  the  network,  whereas,  the  time  increase 
should  be  relatively  small  for  EQFLO  because  the  above  operations  are 
performed  very  efficiently  using  labelling  procedures  on  the  rooted 
spanning  tree. 


Table  4  About  Here 


The  75  side  constraints  made  the  problems  approximately  three  times 
harder.  That  is,  the  pure  networks  were  solved  in  693  seconds  while  it 
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required  1973  seconds  to  solve  the  equal  flow  problem.  Klingman, 
Napier,  and  Stutz  [22]  solved  the  same  35  pure  network  problems  in 
approximately  200  seconds  using  an  advance  start  on  a  CDC  6600  at  the 
University  of  Texas  at  Austin.  Richard  Barr's  best  time  on  these 
35  test  problems  Is  104  seconds  using  ARC  II  [3].  This  difference  in 
time  is  due  to  the  fact  that  EQFLO  is  a  real  code  (as  opposed  to  all- 
integer),  uses  an  all  artificial  start,  and  does  not  use  the  advanced 
data  structure  or  candidate  list  incorporated  in  ARC  II. 

These  35  problems  were  the  largest  that  could  be  solved  using 
NETSIDE  under  a  core  storage  limitation  of  301K  octal  words.  However, 
EQFLO  required  much  less  storage;  approximately  50%  less  than  NETSIDE. 
This  additional  storage  for  NETSIDE  results  from  the  working  basis 
inverse  and  the  arrays  required  during  the  reinversion  process. 
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VII.  SUMMARY  AND  CONCLUSIONS 


This  paper  presents  a  new  procedure  for  the  equal  flow  problem. 
Unlike  the  simplex  method  for  the  network  problem  with  side  constraints, 
this  new  procedure  does  not  require  a  working  basis.  Ve  have  showed 
that  using  the  subgradient  optimization  technique,  the  equal  flow 
problem  may  be  solved  as  two  sequences  of  pure  network  problems.  One 
sequence  corresponds  to  a  lower  bound  while  the  other  corresponds  to 
an  upper  bound.  In  the  lower  bound,  each  network  differs  from  the 
previous  one  in  that  the  cost  vector  has  changed.  In  the  upper  bound, 
each  network  differs  from  the  previous  one  in  that  the  right  hand  side 
has  changed.  While  solving  the  pure  network  problems  with  these 
changes  in  the  problem  data,  a  reoptimization  procedure  is  used  to 
obtain  a  good  starting  solution.  Our  technique  terminates  when  the 
difference  between  two  bounds  is  within  a  prespecified  tolerance. 

Subgradients  for  upper  bounds  are  computed  using  the  optimal 
dual  variables  obtained  by  solving  the  pure  network  problems.  The  sub- 
gradients  for  lower  bounds  are  the  difference  between  the  flows  on  the 
equal  flow  arcs,  obtained  while  solving  the  Lagrangean  relaxation.  The 
projection  operation  is  easily  implemented.  The  step  sizes  (i)  and 
(Hi),  described  in  Section  II,  are  used  for  lower  and  upper  bounds, 
respectively.  For  these  step  sizes,  we  are  guaranteed  a  solution  at 
which  the  objective  value  is  arbitrarily  close  to  the  optimal  objective 
value. 

We  solved  all  test  problems  twice;  once  with  the  same  step  6ize  for 
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all  equal  flow  pairs,  and  once  with  different  step  sizes  for  each  pair. 
The  tests  were  conducted  on  a  set  of  35  randomly  generated  problems  and 
a  comparison  was  made  with  NETSIDE,  a  code  that  is  designed  to  solve 
network  problems  with  side  constraints.  On  the  average,  our  code  ran 
approximately  3  times  faster.  However,  it's  performance  improved 
substantially  as  the  problem  size  increased.  The  new  algorithm  requires 
only  50%  of  the  core  storage  required  by  NETSIDE. 
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Table  1 


NETGEN  Test  Problems 


Problem  Number  Number  of  Nodes 


Transportation  Problems 


1 

100  X 

100 

2 

100  X 

100 

3 

100  X 

100 

4 

100  X 

100 

5 

100  X 

100 

6 

150  X 

150 

7 

150  X 

150 

8 

150  X 

150 

9 

150  X 

150 

10 

150  X 

150 

ent  Problems 

11 

200  X 

200 

12 

200  X 

200 

13 

200  X 

200 

14 

200  X 

200 

15 

200  X 

200 

Capacitated  Network  Problems 


16 

400 

17 

400 

18 

400 

19 

400 

20 

400 

21 

400 

22 

400 

23 

400 

24 

400 

25 

400 

26 

400 

27 

400 

Uncapacitated  Network  Problems 

28 

1000 

29 

1000 

30 

1000 

31 

1000 

32 

1500 

33 

1500 

34 

1500 

35 

1500 

Number  of  Arcs 


1511 

1700 

2207 

2405 

3100 

3450 

4800 

5470 

6395 

6611 


1900 

2650 

3400 

4150 

4900 


1374 

2511 

1374 

2511 

1484 

2904 

1484 

2904 

1398 

2692 

1398 

2692 


3000 

3500 

4500 

4900 

4492 

4535 

5257 

5880 
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Table  2  Comparison  of  NETSIDE  and  EQFLO  on  35  Test  Problems 
(Same  step  size  for  every  equal  flow  pair) 


Problem 

Number 

NETSIDE 

EQFLO 

Optimal  Total 

Objective  Time 

Time 

X  of  Optimal 
at  Termination 

Total 

Pure 

Network 

Lower 

Bound 

Upper 

Bound 

Lower  Upper 
Bound  Bound 

1 

2694547 

30 

50 

7 

19 

24 

98 

108 

2 

2350637 

24 

58 

7 

23 

28 

95 

105 

3 

1939836 

27 

127 

9 

55 

63 

99 

110 

4 

1612265 

33 

79 

10 

33 

36 

98 

108 

5 

1480741 

33 

40 

12 

13 

15 

97 

108 

6 

2472907 

71 

46 

22 

9 

15 

98 

108 

7 

2236784 

96 

59 

28 

14 

17 

97 

107 

8 

2223900 

84 

58 

32 

11 

15 

99 

108 

9 

1839835 

115 

44 

36 

6 

2 

98 

105 

10 

2291942 

105 

79 

36 

19 

24 

96 

106 

11* 

4992 

135 

51 

17 

11 

23 

98 

108 

12 

3573 

105 

93 

23 

20 

50 

95 

105 

13 

3142 

103 

78 

27 

16 

35 

98 

108 

14 

2787 

118 

34 

31 

1 

2 

99 

101 

15 

2795 

150 

127 

35 

31 

61 

97 

108 

16 

82161432 

43 

8 

6 

1 

1 

99 

107 

17 

45601025 

66 

13 

8 

4 

1 

99 

105 

18 

81600312 

40 

8 

6 

1 

1 

99 

106 

19 

45601025 

66 

12 

8 

3 

1 

99 

102 

20 

74065202 

40 

9 

6 

2 

1 

99 

108 

21 

40137087 

44 

11 

8 

2 

1 

99 

101 

22 

73429862 

32 

8 

6 

1 

1 

99 

109 

23 

39354594 

33 

11 

8 

2 

1 

99 

101 

24 

85926653 

91 

7 

3 

3 

1 

98 

104 

25 

58203746 

66 

9 

5 

3 

1 

99 

101 

26 

74267081 

65 

6 

3 

2 

1 

97 

102 

27 

47295659 

57 

7 

4 

2 

1 

99 

107 

28 

131316225 

201 

31 

20 

9 

2 

99 

107 

29 

113594497 

260 

167 

25 

72 

70 

98 

107 

30 

90569484 

337 

243 

23 

111 

109 

91 

106 

31 

84943754 

296 

44 

24 

16 

4 

99 

109 

32 

180390305 

529 

80 

48 

25 

7 

98 

109 

33 

205246112 

453 

83 

47 

23 

13 

98 

108 

34 

166247998 

477 

95 

51 

24 

20 

96 

106 

35 

A  _ 

163964307 

_  1  A 

503 

68 

52 

11 

5 

99 

107 

*  P0  •  10. 
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Table  3  Comparison  of  NETSIDE  and  EQFLO  on  35  Test  Problems 
(Different  step  sizes  for  different  pairs) 


Problem  NETSIDE  EQFLO 

Number 


Optimal 

Objective 

Total 

Time 

Time 

%  of  Optimal 
at  Termination 

Total 

Pure 

Network 

Lower 

Bound 

Upper 

Bound 

Lower 

Bound 

Upper 

Bound 

1 

2694547 

30 

47 

7 

16 

24 

97 

108 

2 

2350637 

24 

50 

7 

19 

24 

94 

105 

3 

1939836 

27 

105 

9 

42 

54 

99 

109 

4 

1612265 

33 

74 

10 

29 

35 

98 

108 

5 

1480741 

33 

37 

12 

10 

15 

95 

106 

6 

2472907 

71 

46 

22 

9 

15 

98 

107 

7 

2236784 

96 

57 

28 

13 

16 

97 

105 

8 

2223900 

84 

42 

32 

4 

6 

98 

109 

9 

1839835 

115 

44 

36 

6 

2 

98 

105 

10 

2291942 

105 

76 

36 

19 

21 

96 

105 

11* 

4992 

135 

53 

17 

9 

27 

94 

104 

12 

3573 

105 

79 

23 

14 

42 

95 

105 

13 

3142 

103 

63 

27 

11 

25 

98 

108 

14 

2787 

118 

34 

31 

1 

2 

99 

101 

15 

2795 

150 

108 

35 

24 

49 

98 

107 

16 

82161432 

43 

8 

6 

1 

1 

99 

107 

17 

45601025 

66 

13 

8 

4 

1 

99 

105 

18 

81600312 

40 

8 

6 

1 

1 

99 

106 

19 

45601025 

66 

12 

8 

3 

1 

99 

102 

20 

74065202 

40 

9 

6 

2 

1 

99 

108 

21 

40137087 

44 

11 

8 

2 

1 

99 

101 

22 

73429862 

32 

8 

6 

1 

1 

99 

109 

23 

39354594 

33 

11 

8 

2 

1 

99 

101 

24 

85926653 

91 

7 

3 

3 

1 

98 

104 

25 

58203746 

66 

9 

5 

3 

1 

99 

101 

26 

74267081 

65 

6 

3 

2 

1 

97 

102 

27 

47295659 

57 

7 

4 

2 

1 

99 

107 

28 

131316225 

201 

31 

20 

9 

2 

99 

107 

29 

113594497 

260 

81 

25 

28 

28 

97 

107 

30 

90569484 

337 

148 

23 

62 

63 

93 

104 

31 

84943754 

296 

44 

24 

16 

4 

99 

109 

32 

180390305 

529 

80 

48 

25 

7 

98 

109 

33 

205246112 

453 

78 

47 

22 

9 

98 

108 

34 

166247998 

477 

86 

51 

23 

12 

97 

107 

35 

163964307 

503 

68 

52 

11 

5 

99 

107 

*  P0  *  10. 
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Table  4  Sunmary  of  Computational  Results 
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1.  THE  SUBGRADIENT  ALGORITHM 


ut  ,  be  .  finite  convex  functional  on  F"'  *«»  *“”• 

define  the  eubdif ferential  of  g  at  y  by: 

3g(y)  -  {n  :  neR*  and  for  all  ztRn,  »W  >  ^(e-y)). 

Any  ncBg(y)  is  called  a  Butadiene  of  g  at  y.  It  *  veil  hnown  that 
„  ,  is  .  point  et  which  g  It  dlf fex.nti.hle.  then  8*(y)  ‘  «’  ' 

singleton  set.  n 

uc  c  p  ,  he  .  closed  end  convex  euheet  of  »  .  **  y<*  ' 

define  the  zrolSSli™  of  V  on  G'  dK“ted  ^  ^ %  ^  ^  ^  ? 

of  6  such  thet  fot  .11  tec.  ||«y)-,ll  <  « ^  " 
the  projection  exists  in  this  ceee  end  thet  fot  .11  x,y«  T  y 

<  H  *-y!l  • 

Coneidet  the  nonlinesr  programing  proble.  given  by. 

(NLP/SD) 

minimize  g(y) 
subject  to  y£G, 

.  *  «n  vc c  3a(v)  1*  *  and  that  the  set  of  optimal 

where  we  assume  that  for  all  y£G,  giy 

points  IM.  we  denote  the  optimal  objective  value  by  1. 

The  eubgt.di.nt  opti.ix.tion  algorithm  for  the  eolution  of  WtP/Sh 
...  fixer  introduced  by  Shot  [111  .nd  may  be  viewed  ..  .  genex.lix.tion 
of  the  eteepeet  de.cent  method  in  which  .ny  eubgx.di.nt  i.  eubetituted  fox 
the  gx.di.nt  .t  .  point  whet,  the  gradient  doe.  not  exiet.  This  .Igorithm 

. . .  of  positive  et.P  .ire.  <»s>,  which  in  turn  depend  on  .  pxe- 

determined  sequence  of  fixed  constent.  U,}  .nd  (in  sum  ce.ee)  c.xt.in 
other  Quantities. 
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SUBGRADIENT  OPTIMIZATION  ALGORITHM 

Step  0  (Initialization) 

Let  yQeG  and  set  i  ■*-  0. . 

Step  1  (Find  Subgradient  and  Step  Size) 

Obtain  some  Ti^eBgtyj). 

If  *  0,  terminate  with  y^  optimal;  otherwise,  select  a  step 
size  . 

Step  2  (Move  to  New  Point) 

Set  yi+1  P(y^-s.r)  ) ,  i  ■*-  i+1,  and  return  to  step  1. 

Unfortunately,  the  termination  criterion  in  step  1  may  not  hold 
at  any  member  of  T  and  is  thus  computationally  ineffective.  Hence,  some 
other  stopping  rule  must  be  devised.  In  practice  this  is  often  a  limit 
on  the  number  of  iterations.  The  functional  values  produced  by  the 
algorithm  will  be  denoted  by  g^  «  g(y^). 

Various  proposals  have  been  offered  for  the  selection  of  the 
step  sizes.  Four  general  schema  which  have  been  suggested  are: 


xi  • 

(i) 

\  /  II n4  II , 

(2) 

\  /  IlnJI 2  , 

(3) 

\  (gt-P)  /  IlnJI 2  , 

(A) 

where  p,  the  target  value,  is  an  estimate  of  y  and  all  0. 

The  papers  of  Polyak  [9]  and  Held,  Wolfe,  and  Crowder  [6]  have 
provided  the  major  impetus  for  widespread  practical  application  of  the 
algorithm.  Schema  (4)  has  proven  to  be  a  particularly  popular  choice 


2 


among  experimenters.  Theorem  4  of  Polyak  [9]  Is  the  most  often  quoted 
convergence  result  Justifying  use  of  this  schema.  For  many  mathematical 
programming  models,  the  target  value  is  a  lower  bound  on  the  optimum 
(i.e.,  [1,  2,  3,  4,  7,  8,  10]).  For  this  case  Polyak’s  Theorem  4,  using 
schema  (4),  requires  that  A^  *  1  for  all  1.  For  all  the  above  studies, 
a  decreasing  sequence  of  A’s  was  found  to  work  better  than  A^  ■  1  for 
all  1.  Hence,  the  existing  theory  did  not  justify  what  we  had  found  to 
work  well  in  practice.  The  objective  of  this  paper  is  to  present  new 
theoretical  results  which  help  to  explain  what  has  been  found  to  work 
well  in  practice.  Specifically,  we  have  generalized  Polyak’s  result  for 
a  decreasing  sequence  of  A's.  In  addition,  we  also  loosen  slightly  the 
restrictions  imposed  on  the  sequence  {A^}  when  the  target  value  is  larger 
than  Y- 
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II.  POLYAK'S  CONVERGENCE  RESULTS 


Some  of  Che  convergence  results  for  Che  subgradient  optimization 
algorithm  appear  unusual  in  that  they  specify  only  that  a  functional 
value  vlthin  a  given  tolerance  of  the  optimal  value  y  will  eventually  be 
produced.  The  results  of  Theorem  4  of  Polyak  [9]  use  the  following 
general  restrictions  on  the  sequence  fx  used  with  schema  (4): 

0  <  a  <  <  B  <  2,  (5) 

where  a  and  £  are  fixed  constants. 

The  results  contained  in  this  theorem  include: 
under  (4),  (5),  and  (essentially)  the  assumption  that  there  is  some 
K>0  such  that  ||  ||  K 

(A)  if  p  >  y,  either 

(a)  there  is  some  n  such  that  g^  <  p, 
or 

(b)  all  gR  >  p  and  lim  g^  *  p; 

and 


(B)  if  p  <  Y  and  all  Xn  ■  1,  given  6  >  0, 

there  is  some  n  such  that  gn  <  y  +  (Y"p)  +  $. 

In  the  next  section  we  will  relax  condition  (5)  to  the  following: 

0  <  -  8  <  2  and  I  X£  -  <*>,  (6) 

where  £  is  a  fixed  constant,  and  we  will  present  a  generalization  of  (B) 
for  a  decreasing  sequence  {X^}. 
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III.  NEW  CONVERGENCE  RESULTS 


The  main  results  In  this  section  appear  in  Propositions  5  and  7. 
Proposition  5  corresponds  to  part  A  of  Polyak's  Theorem  4  with  slightly 
weaker  conditions  on  the  sequence  {X^},  and  Proposition  7  is  a  generali¬ 
zation  of  part  B  of  Theorem  4. 

Proposition  1 

If  yeG,  then 

II  y-yi+1  II 2  <  II  y-yt  II 2  +  s2  ||  n1ll2  +  2si  (g(y)-gi). 

Proof 

Let  yeG. 

II  y-yi+i  II 2  *  II  y-P(yi-sin1)  ||2 

-  II  P(y)>P(y1-s1n1)||  2 

<  ||  y-yi+sjLrii||  2 

-  II  y-y± II  2+s2||  rtjl  2+2s1ni  •  (y-y^ 

<  II  y-y±  1 1  2+s2||  njl  2+2s1(g(y)-g1). 

Proposition  2 

If  yer,  then  under  (4), 

lly-yi+i  II 2  i\\  y^ll  (gi-p)l>'i(gi-p)-2(gi-Y)]/||  rii|| 2. 

Proof 

Let  yer.  Substituting  in  Proposition  1  for  Sj^  from  (4)  and 
using  g(y)*Y,  we  obtain 

II  y-yi+i  II 2  f  II y-yA  ll2+*2  (gi-p)2/lhill  2+2XA  (g1-p)  (Y-g^  /  1!  n±ll 2 
■  II y-y±  I! 2+xi  (gj-p)!^  (gj-p^Cg^Y)]  /  II  ni ||  2. 
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Proposition  3 

If  yeP,  p>Y,  and  >  p,  then  under  (4), 

II  y_yi+lH  2  -  Hy“yiH  2  +  Xi  (Xi"2>  <8i“P>2  /II  ^11  2‘ 

Proof 

Let  yer,  p>y,  and  g±>p  . 

Now,  \±  (gi-p)-2(gi-y)  <  (g±-p)-2(g±-p)  -  (A±-2>  (g^p). 

Thus,  x±  (gi-p)[x1(g1-p)-2(gjL-Y)]  /  IlnJI2  <  x±(Xr2)  (gl-p)2/ IlnJI  2. 

The  result  now  follows  from  Proposition  2. 

Proposition  4 

If  yer,  p>Y,  and  all  g^p,  then  under  (4)  with  all  X^2,  there 

2 

is  some  ^  such  that  lim  Hy-y^H  *  ip. 

Proof 

Let  yer,  p>Y»  all  X^<2,  and  all  g^>p.  Since  each  X^<2,  then 

2  2 
also  each  X±  <X±— 2)  (g^p)  /||  njl  <  0,  and  from  Proposition  3,  {  ||  y-yjl  ) 

is  a  monotone  nonincreasing  sequence.  This  sequence  is  bounded  below  by 

zero  and  thus  converges  to  some  value,  say  \J>. 

Proposition  5 

If  p>Y  and  there  is  some  K>0  such  that  all  ||  ri^||  <tc,  then  under 
(4)  and  (6),  given  6>0,  there  is  some  M  such  that  g^j  <  p+d. 

Proof 

Let  6>0  be  given,  with  p>Y*  and  all  ||  t^H  «•  Suppose,  contrary 
to  the  desired  result,  that  all  gj^  >  p+6.  Take  any  yer.  Then  from 
Proposition  3, 

xi(2-x1)(gi-p)2/||  n1ll  2  <  ll  y-y1  il  2-l!  y-y±+1  II  2. 
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(7) 


Since  A  <  8  <  2,  ||  n± ||  «,  and  gi~  p  >  6, 

\  (2-B)S2/<2  <  ||  y-yjl  2  -  ||y-yi+1!l  2. 

Adding  together  the  inequalities  obtained  from  (7)  by  letting  i  take  on 
all  values  from  0  to  n,  we  obtain 

(V*  *  ,+xn)  tf-w*2/*2  ^  h-y0  II 2  -  lly-yn+i"2*  (8) 

As  n  goes  to  %  the  left  side  of  (8)  goes  to  «°,  whereas,  by  Proposition  A, 

2  2 

the  right  side  of  (8)  goes  to  ||y-y0  ||  »  a  contradiction. 

Proposition  5  gives  a  practical  convergence  result  when  the  target  exceeds 
the  optimal  value.  At  worst  we  eventually  obtain  an  objective  value 
arbitrarily  close  to  the  target  value. 

Proposition  6 

If  yer,  g^P,  and  .<  6/2,  then  under  (A), 

lly-y1+ill 2  i  Ily-yJI 2 

+  \(g1-P)(2-6)  [  (Y-gi)  +  (8/(2-g))  (y-p)  3  /  lln.ll2. 

Proof 

Let  yer,  g^P,  and  Ai  <  8/2. 

Now,  (gi-p)-2(g1-y)  <  B(gi-p)-2(gi-y) 

«  B(gi-p)-(2-B) (gi-y)-B(g1-y) 

•  6(y-p)+(2-6) (y-g±) 

-  (2-B)  l  (Y-g^  +  (8/ (2-B))  (Y-P)  3. 

Thus,  Xi  (gi-p)  [Xi(gi-p)-2(gi-Y)  3  /  H  nj)  2 

1  \  (g^P)  (2-B)  l  (Y-gj)  +  (8/ (2-8) )  (Y-P)  3  /  linjl2. 

The  result  now  follows  from  Proposition  2. 

Proposition  7 

If  p<Y  and  there  is  some  tc>0  such  that  all  ||  n^||  «,  then  under 
(A)  and  (6),  given  6>0,  there  is  some  M  such  that  <  y  +  (6/ (2— B) )  (y-p)  + 
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Proof 


Let  5>0  be  given,  with  p<Y»  and  all  ll  n±||  «.  Suppose,  contrary 
to  the  desired  result,  that  all  >  Y+(6/(2-6))  (y~p)  +  d,  or 
(Y-g^)  +  (6/(2-B))  (y~p)  <  "*•  Since  B<2  and  g^p,  then 

\  (g  1-P)  (2-0)  [  (Y-gi)  +  (B/(2-B))  (Y-P)  1  /  ||n1If2 
<  -6X1  (8l-p)  (2-0)  /  ||  niH  2.  (9) 

Take  any  yer.  Then  by  (9)  and  Proposition  6,  we  have  that 

6Ai(gi"p)  (2_e)  /  'I  niH  2  <  II  y“yiH  2“H  y“yi+lH  2‘ 

Since  ||  n .  H  <<  and  gi>Y>P»  then  also 

X^Y-P)  (2-6)  /  K2  <  II  y-y±||  2-||  y-yi+1H  2.  (10) 

Adding  together  the  inequalities  obtained  from  (10)  by  letting  i  take  on  all 
values  from  0  to  n,  we  obtain 

( V  *  *+V  6  (Y"P)  (2~6)  '  <2  <  II  y'y0  II  2“il  y-yn+l  II  2-  (11) 

As  n  goes  to  °°,  the  left  side  of  (11)  goes  to  °°,  whereas,  by  Proposition  4, 

2  2 

the  right  side  of  (11)  goes  to  ||  y-y ^  ||  -ip  ,  a  contradiction. 

The  above  is  our  generalization  of  Polyak’s  Theorem  4  Part  B. 

At  worst  we  eventually  obtain  an  objective  value  whose  error  is  arbitrarily 
close  to  0/ (2-6)  times  the  error  present  in  the  target  value  estimate  of  y. 


IV.  CONCLUSIONS 


Proposition  5  gives  the  convergence  result  obtained  under  (4) 
and  (6)  for  a  target  value  at  or  above  the  optimal  value.  It  is 
readily  apparent  that  Proposition  5  is  compatible  with  Polyak's  result 
(A).  Proposition  7  gives  the  corresponding  result  for  a  target  value 
under  the  optimal  value.  We  have  found  this  to  be  a  more  practical 
result  (see  e.g.,  (1,  2,  3,  4,  7,  8,  10]).  Taking  B»l,  we  have 
Polyak's  result  (B)  as  a  special  case  of  Proposition  7.  Proposition  7 
shows  more  clearly  the  dependence  of  the  demonstrably  attainable  error 
on  the  upper  bound  8  for  {X^}.  This  paper  has  not  addressed  the  question 
of  any  convergence  rate  associated  with  the  use  of  (4)  and  (6).  Goffin 
15]  has  provided  such  results  when  schema  (2)  is  used. 
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ABSTRACT 

An  important  class  of  mathematical  programming  models  which 
are  frequently  used  in  logistics  studies  is  the  model  of  a  network 
problem  having  additional  linear  constraints.  A  specialization  of 
the  primal  simplex  algorithm  which  exploits  the  network  structure 
can  be  applied  to  this  problem  class.  This  specialization  maintains 
the  basis  as  a  rooted  spanning  tree  and  a  general  matrix  called  the 
working  basis.  This  paper  presents  the  algorithms  which  may  be  used 
to  maintain  the  inverse  of  this  working  basis  as  an  LU  factorization, 
which  is  the  industry  standard  for  general  linear  programming  soft¬ 
ware.  Our  specialized  code  exploits  not  only  the  network  structure 
but  also  the  sparsity  characteristics  of  the  working  basis.  Compu¬ 
tational  experimentation  indicates  that  our  LU  implementation  results 
in  a  50%  savings  in  the  nonzero  elements  in  the  eta  file,  and  our 
computer  codes  are  approximately  twice  as  fast  as  MINOS  and  XMP  on  a 
set  of  randomly  generated  multicommodity  network  flow  problems. 
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I.  INTRODUCTION 


Good  software  for  solving  linear  programming  models  is  one  of  the 
most  important  tools  available  to  the  logistics  engineer.  For  logistics 
studies,  these  linear  programs  frequently  involve  a  very  large  network 
of  nodes  and  arcs,  which  may  be  duplicated  by  time  period.  For  example, 
nodes  may  represent  given  cities  at  a  particular  point  in  time  while 
arcs  represent  roads,  railways,  and  legs  of  flights  connecting  these 
cities.  Some  nodes  are  designated  as  supply  nodes,  others  demand  nodes, 
while  some  may  simply  represent  points  of  transshipment.  The  mathematical 
model  characterizes  a  solution  such  that  the  supply  is  shipped  to  the 
demand  nodes  at  least  cost  while  not  violating  either  the  upper  or  lower 
bounds  on  the  flow  over  an  arc. 

If  the  main  structure  of  a  logisitcs  problem  can  be  captured  in  a 
network  model,  then  the  size  of  solvable  problems  becomes  enormous. 

Hence,  more  realistic  situations  can  be  modelled  that  would  otherwise  lie 
outside  the  domain  of  general  linear  programming  techniques.  For  example, 
one  current  logistics  planning  model  involves  200  nodes  and  (365  days/yr) 
(30  years)  *  10,950  time  periods  to  give  over  2,000,000  constraints. 
Network  problems  having  20,000  constraints  and  20,000,000  variables  are 
solved  routinely  at  the  U.  S.  Treasury  Department. 

Unfortunately,  the  pure  network  structure  may  require  simplification 
of  the  problem  to  the  point  that  key  policy  restrictions  must  be  omitted. 
The  work  presented  in  this  study  builds  upon  existing  large-scale  network 
solution  technology  to  allow  for  the  inclusion  of  arbitrary  additional 


constraints.  Typical  constraints  include  capacities  on  vehicles 
carrying  different  types  of  goods,  restrictions  on  the  total  number 
of  vehicles  available  for  assignment,  and  budget  restrictions.  The 
addition  of  even  a  few  non-network  constraints  can  greatly  enhance 
the  realism  and  usability  of  these  models.  Our  approach  exploits — 
to  as  great  an  extent  as  possible — the  traditional  network  portion  of 
the  problem  while  simultaneously  enforcing  any  additional  restrictions 
imposed  by  the  practitioner. 

For  general  linear  programming  systems,  the  most  important 
component  is  the  algorithm  used  to  update  the  basis  inverse.  Due 
to  the  excellent  sparcity  and  numerical  stability  characteristics, 
an  LU  factorization  with  either  a  Bartels-Golub  or  Forrest-Tomlin 
update  has  been  adopted  for  modern  linear  programming  systems.  For 
pure  network  problems,  the  basis  is  always  triangular  and  corresponds 
to  a  rooted  spanning  tree.  The  modern  network  codes  which  exploit 
this  structure  have  been  found  to  be  from  one  to  two  orders  of 
magnitude  faster  than  the  general  linear  programming  systems.  In 
this  paper,  we  have  combined  these  two  powerful  techniques  into  an 
algorithm  for  solving  network  models  having  additional  side  constraints. 

Let  A  be  an  m  x  n  matrix,  let  c  and  u  be  n-component  vectors,  and 
let  b  be  an  m-component  vector.  Without  loss  of  generality,  the  linear 
program  may  be  stated  mathematically  as  follows: 

minimize  c  x  (1) 
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subject  to: 


A  x  *=  b 


(2) 


0  <  x  <  u  .  (3) 

The  network  with  side  constraint  model  is  a  special  case  of  (1)  -  (3) 
in  which  A  takes  the  form 


A  - 

M  i 

I 

S  I  P 

* 

—  — 

-  , 

where  M  is  a  node-arc  incidence  matrix. 


1.1  Applications 

There  are  numerous  applications  of  the  network  with  side  constraint 
model.  Professor  Glover  and  his  colleagues  have  solved  a  large  passenger- 
mix  model  for  Frontier  Airlines  and  a  large  land  management  model  for 
the  Bureau  of  Land  Management  (see  [  7,  8]).  A  world  grain  export  model 
has  been  solved  to  help  analyze  the  port  capacity  of  U.  S.  ports  during 
the  next  decade  (see  [2]).  A  cargo  routing  model  is  being  used  by  the 
Air  Force  Logistics  Command  to  assist  in  routing  cargo  planes  for  the 

distribution  of  serviceable  spares  (see  [l]).  Lt.  Col.  Dennis  McLain, 
has  developed  a  large  model  to  assist  in  the  development  of  a 
casualty  evacuation  plan  in  the  event  of  a  European  conflict  (see 
[14]).  A  National  Forest  Management  Model  has  been  developed  to 
aid  forest  managers  in  long  term  planning  for  national  forests 


(see  [10]).  In  addition,  work  is  currently  underway  which  attempts  to 
convert  general  linear  programs  into  the  network  with  side  constraint 
model  (see  [4,  16]). 

1.2  Objective  of  In\ ^stigation 

Due  to  both  storage  and  time  considerations,  the  basis  inverse  is 
maintained  as  an  LU  factorization  in  modern  LP  software  (see  [3,5,  15]). 
The  objective  of  this  investigation  is  to  extend  these  ideas  to  the 
primal  partitioning  algorithm  when  applied  to  the  network  with  side 
constraints  model. 


1.3  Notation 

The  ith  component  of  the  vector  a  will  be  denoted  by  a^.  The 

(i,j)th  element  of  the  matrix  A  is  denoted  by  A„.  A(i)  and  A[i] 

denotes  the  it^1  column  and  i^  row  of  the  matrix  A,  respectively. 

v 

0  denotes  a  vector  of  zeroes,  1  denotes  a  vector  of  ones,  and  e  denotes 
a  vector  with  a  1  in  the  k*1*1  position  and  zeroes  elsewhere.  Sigma  is 
used  to  denote  the  scalar  signum  function  defined  by 


1,  if  y  >  0 

o(y)  *  /  0,  if  y  ■  0 

-1,  if  y  <  0  . 
V. 

The  identity  matrix  is  given  by  "I". 
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II.  THE  PRIMAL  SIMPLEX  ALGORITHM 


We  assume  that  A  has  full  row  rank  and  that  there  exist  a  feasible 
solution  for  (1)  -  (3).  Given  a  basic  feasible  solution,  we  may  par¬ 
tition  A,  c,  x,  and  u  into  basic  and  nonbasic  components,  that  is, 

A  =  [B  !  N],  c  =  [cB  !  cN],  x  ■  [xB  j  x^],  and  u  *  [uB  j  uN] .  Using  the 
above  partitioning,  the  primal  simplex  algorithm  may  be  stated  as 
follows : 


PRIMAL  SIMPLEX  ALGORITHM 


0.  Initialization.  Let  [xB  !  xN]  be  a  basic  feasible  solution. 

TJ  _-1 

1.  Pricing.  Let  II  *  c  B  .  Define 

*=  {i  :  *  0  and  II  N(i)  >  c^}, 

$2  m  (i  :  x^  =  u^  and  II  N(i)  <  c^}. 

If  iJ^U  ^  >  terminate  with  [xB  j  xN]  optimal;  otherwise, 

select  k  e  ^  (j  and  set  5  +  1  if  k  £  ^  and  6  -1,  otherwise. 

2.  Ratio  Tza. C.  Set  i  *■  B-1  N(k).  Set 

B 


min 


A1  -  a<yJ)-i  (  ,y< 


J 


min 

A2  *►  -o(y^)=6 


B  B 

u ,  -  X. 

J  3 


N 

Set  A  min  {A^,  u^}, 


If  A  y  °°,  then  go  to  3;  otherwise,  terminate  with  the  conclusion 
that  the  problem  is  unbounded. 
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Select  any  1  e  y  In  the  basis,  replace  B (2)  with  N(k), 

update  the  inverse  of  the  new  basis,  and  return  to  step  1, 
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III.  THE  PARTITIONED  BASIS 


The  network  with  side  constraint  model  may  be  stated  as 
follows : 


minimize 

11,22 
c  x  +  c  X 

(4) 

subject  to: 

M  xl  * 

b1 

(5) 

S  x1  +  P  x2  *= 

b2 

(6) 

0  <  x1  <  u1 

(7) 

2  2 

0  <  xZ  <  u  , 

• 

(8) 

We  may  assume  without  loss  of  generality  that, 

(i)  The  graph  associated  with  M  has  n  nodes  and  is  connected  (i.e. 
there  exists  an  undirected  path  between  every  pair  of  nodes). 

(ii)  [S  ;  P]  has  full  row  rank  (i.e.,  rank  [S  j  P]  =  m) . 

(iii)  Total  supply  equals  total  demand  (i.e.,  1  =  0). 

Since  the  rank- of  system  (5)  is  one  less  than  the  number  of 
rows,  we  add  what  has  been  called  the  root  arc  to  (5)  to  obtain 

M  x"*"  +  e15  a  «  b^" 

where  0  <  a  <  0  and  1  <  p  <  n. 

Then  the  constraint  matrix  for  the  network  with  side  constraints  model 
becomes 


It  is  well-known  that  every  basis  for  A  may  be  placed  in  the 


form 


B  = 


where  T  corresponds  to  a  rooted  spanning  tree  and 


(9) 


T  1  +  T_1  C  Q'1  D  T_1 


-T-1  C  Q_1 


-Q-1  D  T*1  }  Q_1 

where  Q«F-DT^C.  The  objective  of  this  paper  is  to  give 
algorithms  which  maintain  Q-1  as  an  LU  factorization. 


(10) 


IV.  THE  INVERSE  UPDATE 


Recall  that  the  partitioned  basis  takes  the  form 
key  nonkey 


Let 


.-1 


-T-1C 


and  let 

B  *  B  L  = 

The  inverse  update  requires  a  technique  for  obtaining  a  new  Q  J'  after 
a  basis  exchange.  Let  B^,  L^  B^  and  Q  denote  the  above  matrices  at 
iteration  i.  Then  we  want  an  expression  for  in  terms  of  . 


The  transformation  takes  the  form 


■ E  Bi 


-l 


(n) 


where  E  is  either  an  elementary  column  matrix  or  a  permutation  matrix. 


Let  E  be  partitioned  to  be  compatible  with  B.  That  is. 


r-i 


n  m 

By  examining  the  (2,2)  partition  of 

-1 


(E4  -  e3  t-1  C)  q:1 


we  obtain 


(12) 
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In  determing  the  updating  formulae,  we  must  examine  two  major 


cases  with  subcases. 

Case  1.  The  leaving  column  is  nonkey.  For  this  case,  E  takes  the  form 


and  (12)  reduces  to  Q^'  “  E^  Q^. 


Case  2.  The  leaving  column  is  key. 


Let  y  *  a**  T  ^  C.  If  #  0,  then  the  k*"*1  column  of  C  can  be  interchanged 


.th 


with  the  j  column  of  T  and  the  new  T  will  be  nonsingular. 
Subca&e.  2a.  y  t  0.  Suppose  Yk  +  0. 


Then  E.  -  T  1  C  reduces  to 
A  3 


'  - 

r- 

I 

i  i 

i  i 

i  i 

1  1 

: 

R  * 

-ej  T_1  C 

row  j 

•  (13) 

- 

i  i 

<  1 

1  I.1  J 

,  and 

<>lii  - R  “l1- 

Case 

1  is 

applied  to  complete  the  update. 

Subca&z  2b.  1 

-  0. 

For 

this  case  no  interchange  is  possible, 

the  entering 

column  becomes  key,  and  Q  ^  =  Q^. 
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1  s 

where  F  ,  . F  are  a  combination  of  row  and  column  etas.  The 
right  side  of  (14)  is  referred  to  as  the  eta  file  where  only  the 
non-identity  rows  and  columns  are  stored.  Suppose  that  the 
column  of  Q  is  replaced  by  §(k)  to  form  the  new  m  by  m  working  basis 
Q.  This  section  presents  algorithms  which  may  be  used  to  update  (14) 

a-1 

to  produce  Q  in  the  same  form* 


5.1  Nonkey  Column  Leaves  The  Basis 


We  will  show  that  Q-1  =  U1  . . .  Um-1  u“  L®  FS  ...  F1. 

V 

If  k  <  m,  then  let  R  =  I  and 

Q-1  ■=  U1  ...  Uk  Rk  Uk+1  ...  Um  FS  ...  F1. 

k.  k+1 

We  next  define  a  new  upper  eta,  0  ,  and  a  new  row  eta,  R  , 


that 


Rk  Uk+1  =  ^  Rk+1 


Substituting  (16)  into  (15)  yields 


-1  TI1  ITk  k+1  TTk+2 
Q  “  U  ...  U  U  R  U 


0°  Fs  ...  F1. 


'Mc+l  k+2 

We  again  define  two  new  eta’s,  U  and  R  ,  such  that 

Rk+1  yk+2  .  £k+l  Rk+2> 

-1  k  'kk  'kk+T 

Substituting  (IS)  into  (17)  yields  Q  «  U  ...  U  U  U  ] 
U™  FS  ...  F1. 


Repeating  this  process  eventually  yields 


n-l  1  k  CVnt— 1  m  s  1 

Q  U  •  •  •  U  U  •  #  •  U  R  F  •  •  •  F  « 


Let  l  -  R"  fs  ...  F1  Q(k),  let 


L» 


-1  T  — 


(15) 

such 

(16) 

(17) 

(18) 

,k+2  k+3 

w  u  •  • 

(19) 


’  -  '••****' *«****► 


I 


and  let 


u"1 


1 

1  -Y 

!  yi 

1  1  •• 

i  -"^k-l 

1  1 

...  j 

1 

1  1 
i  _ 

i 

i 

i 

i 

i 

i 

i 

i 

i 

i 

i 

Then  Um  L®  J 
>RffiFS  ... 


k  ^—1 

e  and  we  will  show  that  Q 


U1  ...  uk_1  uk 


U1" 


We  now  present  the  algorithm  which  updates  the  LU  representation 
of  Q  1  when  the  leaving  column  is  nonkey.  Assume  that  Q(k)  is  replacing 
Q(k)  in  the  working  basis. 
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A LG  1 :  LU  UPDATE  FOR  MONKEY  LEAVING  COLUMN 


Set  6  ♦  FS  ...  F1  Q(k) . 


If  k  i*  m,  set  £  k,  R  ^  I,  go  to  4. 


Set  Im  I,  where  I  is  m  by  m. 


Set  L®  -  1/6  . 

mm  m 


Set  0m  I,  where  I  is  m  by  m. 

Set  tf0  •* — 6.*  for  1  <  j  <  m. 

J®  j’ 

Stop  with  Q-1  -=  U1  ...  U®-1  if1  Em  FS 
Set  or«-R£[k]  U£+1(£+l). 

Set  R^+1  «-  R£. 

Set  CL  *  “• 

Set  U i  *■  V £+1. 


Set  \,£+l  -  °* 

(R£  U£+1  *  0 Z  R£+1) 


Set  £■*-.£  +  1. 

If  £  <  m,  go  to  4. 

<uk+1  ...  ii”-  8k  ...  tf1  *".) 

Set  6  ♦  RmB. 

Set  Zm  ■*-  I,  where  I  is  m  by  m. 
Set  ^k  ^  1/Bk* 

Set  L1?,  -B./6,  ,  for  k  <  j  <  m. 

jk  J  k’  - 

Set  T5m  I,  where  I  is  m  by  m. 

Set  tJ1?,  -6  ,  for  1  <  j  <  k. 

jk.  j  -  - 


Stop  with  Q-1  -  U1  ...  Uk_1  0k  t5k+1 


D“  Zm  Rm  FS  ...  F1. 


We  now  present  the  justification  for  step  3  of  ALG  1.  For  k  *  m, 
we  claim  that  0  1  *=  U1  ...  U™-1  0m  Fs  F1.  Note  that  Q-1  Q(m)  - 

. . .  Um  ^  tJm  Im  8.  But  by  construction  0™  Em  6  ■  em.  Consider 
Proposition  1. 

Let  8  be  any  m-vector  and  be  any  column  eta.  If  8^  “  0,  then 
E1  6  «  8.  _ _ _ 

1  m-1  m  m  a— 1  a  jj 

By  Proposition  1,  U  . ..  U  e  *  e  .  Therefore,  Q  Q(m)  =  e  .  For 

s  1 

1  <  £  <  m,  let  Y  «  F  ...  F  Q(£).  By  construction  Yj  “  0  for  £  <  j  <  m 

and  Yj£  **  1»  By  Proposition  1,  u*+1  ...  iT1  0“  tm  x  m  y.  By  the  con- 

1  o  1  £  £ 

struction  of  U  ...  IT,  we  have  U  ...U  X  “  -  '  Therefore,  if  the 

leaving  column  is  Q(m),  then  step  3  of  ALG  1  produces  Q 

We  now  present  a  theoretical  justification  for  step  4  of  ALG  1. 

Proposition  2. 


and  R 


row  £ 


column  £ 


* 

where  £  J  £  . 
If 


+ 

column  £ 


-16- 


Hence,  ALG  1  produces  the  updated  working  basis  inverse. 


T 


5.2  Key  Column  Leaves  The  Basis 


In  this  section,  we  present  an  algorithm  for  updating  the  working 
basis  inverse  to  accomplish  a  switch  between  a  key  column  and  a  nonkey 
column.  That  is,  0  =  R  Q  ^  where  R  is  given  by  (IB''  and 

Q-1  =  U1  ...  Um  FS  ...  F1.  (20) 


/V— 1 

We  wish  to  obtain  Q  in  the  same  form  as  (20). 

To  accomplish  this  update,  we  begin  with  Q  ^  =  RU^  ...  U111  FS  ...  F^. 

1  ~-i  .,12  2 

We  apply  Proposition  2  to  R  U  creating  the  factorization  Q  =  U  R  U  ... 

u“  Fs  ...  F1.  We  continue  with  the  application  of  Proposition  2  until 

we  obtain  Q  ^  «  0^  ...  0^  ^  R^  ...  U™  FS  ...  F^.  Proposition  2  does  not 

apply  to  R  U  .  However,  a  simple  update  would  be  to  let  tj  =...=  0®  =  I 

and  use  the  below  factorization: 


2>~1  ?%1  mIQ  „k  k  m  s  _1 

Q  sU  •  • .  U  R  U  . . .  U  F  . • ■  F  . 


LEFT  FILE 


RIGHT  FILE 


This  update  simply  involves  application  of  Proposition  2  until  it  does 
not  apply  (£  -  £*)  and  then  shifting  the  remainder  of  the  left  file 
to  the  right  file.  We  call  this  update  the  TYPE  1  UPDATE. 

We  will  now  give  an  update  in  which  R^  ...  Um  is  modified  as 
opposed  to  moving  them  to  the  right  file.  Let 


* 


-1R- 


■*-  row  k 


,  , .  .  «.k+l  ,  „k+l  .  ,  _k  k+1  «k+l  _k+l 

Then  we  define  matrices  U  and  E  such  that  E  U  =  U  E 


k+1  p,m  _nri-l 


Following  this  procedure,  R  U  ...  U™  can  be  replaced  by  0  ...ft  E 


so  that 


Q-1  -  61  ...  0k_1  ok+1  ...  0"  E”"1  Fs  ...  F1. 

Further,  we  define  a  row  eta  8.  and  a  column  eta  £  such  that  =  ft  ?. 

Therefore, 

a-1  -  o1  ...  ok"1  ok+1  ...  0m  8  P  Fs  ...  F1  . 


left  file 


RIGHT  FILE 


We  call  this  update  the  TYPE  1  UPDATE. 

We  now  present  a  set  of  propositions  which  justify  the  TYPE  2 

UPDATE.  _  _  _  .  _ _ 

Proposition  7. 


V_!  _Yji*r--Yn 
~mi  1 


where  Z  ¥  Z*  and  U 


where 


Y  n,  if  i  *  £ 


The  following  proposition  is  used  to  replace  the  cross  matrix 
m+1  *-■ 

E  with  a  row  eta  R  and  a  column  eta  R.  _____ 

Proposition  8  . 


gijjjgttB 

h  ! 

;  1  » 

vi  | 

V*-y£-i  j 

Y£  !  Y£+l*’*Yn 

1 

l 

l 

l 

y£+l  i 

0  i 

:  I  i 

#  • 

l 

I 

y_  j 

l 

1 

n  i 

L 

Y1'**Y£-1  {  X 


'£+1*  *  *  n 


! 

1 

1 

i  ! 

1 

,  ,T  _ 

yi 

• 

• 

y£-l 

Y 

- 

0 

y£+l 

yn 

where  X  and  Y  are  such  that  XY  =  Yp  “  l  Y.y.> 

^  11 


then  E  *  R  F. 

We  now  present  the  update  algorithm  for  the  case  in  which  the 
column  of  T  is  being  switched  with  the  column  of  C. 


V 


A LG  2:  LU  UPDATE  FOR  A  KEY  LEADING  COLUMN 

1.  Set  R1  ♦  I. 

Set  R1[k]  «-  Y. 

Set  i  ■+■  1. 


2.  If  i  *  k,  go  to  4. 
Set  a  ^  Ri[k]  U^i). 
Set  Ri+1  «-  R1. 

Set  R^1  •»-  a. 

Set  O1  U1. 

Set  0.1.  «-  0. 
kx 


3. 


Set  i  *■  i  +  1  and  go  to  2. 


4.  Set  0 


I. 


„  _k  _  _k  k 
Set  E  R  U  . 


5.  Apply  Proposition  7  to  E^  to  form 

Set  i  i  +  1. 

6.  If  i  <  m,  go  to  5. 

7.  Apply  Proposition  g  to  Em  to  obtain  R  ?  where  X  =  1. 

At  the  completion  of  step  7  we  have  Q  ^  *  0^  ...  R  R  FS 6 7  ...  f\ 
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VI.  COMPUTATIONAL  EXPERIMENTATION 


Three  test  problems  were  selected  for  the  experiment. 

SC205  is  a  staircase  linear  program  which  was  generated  by  Ho  and 
Loute  [  12]  and  transformed  into  a  network  with  side  constraints t 
Gif ford-Pinchot  is  a  model  of  the  Gif ford-Pinchot  National  Forest  [10] 
which  has  also  been  transformed  into  a  network  with  side  constraints. 

RAN  is  a  randomly  generated  problem. 

These  problems  were  first  solved  and  the  pivot  agenda  was  saved. 

That  is,  entering  and  leaving  columns  for  each  pivot  were  saved  on  a  file. 
This  file  was  then  used  by  each  code  so  that  all  three  basis  updates 
follow  the  same  path  to  the  optimum.  The  number  of  nonzeroes  re¬ 
quired  to  represent  Q  *  at  various  points  in  the  solution  process  is 
illustrated  in  Figures  1  and  2.  For  both  problems,  the  LU  Type  2  up¬ 
date  dominated  both  the  LU  Type  1  update  and  the  product-form  code  in 
terms  of  ncnzeroes  in  the  inverse.  The  average  core  storage  required 
for  Q  ^  using  the  product-form  update  is  approximately  double  that 
required  for  the  best  LU  update. 


Figures  1  and  2  About  Here 

Given  the  above  results,  we  developed  three  specialized  network 
with  side  constraints  codes  and  computationally  compared  them  with  three 
general  in-core  LP  systems  and  a  special  system  for  multicommodity  network 
flow  problems.  All  codes  are  written  in  FORTRAN  and  have  not  been 
tailored  to  either  our  equipment  or  our  FORTRAN  compiler.  None  of 

the  codes  were  tuned  for  our  problem  set.  A  brief  description  of  each 
code  follows. 
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NETSIDE1,  NETSIDE2  AND  NETSIDE3  are  our  specialized  network  with 
side  constraints  systems.  The  first  maintains  Q-1  in  product  form,  while 
the  second  and  third  maintain  Q  ^  in  LU  form  using  a  Type  1  and  Type  2 
update,  respectively.  All  use  the  Hellerman  and  Rarick  [11]  reinversion 
routine.  The  working  basis  is  reinverted  every  60  iterations.  The 
pricing  routine  uses  a  candidate  list  of  size  6  with  block  size  of  200. 

MINOS  [15]  stands  for  "a  Modular  In-Core  Nonlinear  Optimization 
System"  and  is  designed  to  solve  problems  of  the  following  form: 

minimize  f(x)  +  cx 

subject  to:  Ax  =  b 

I  <  x  <  u 

where  f(x)  is  continuously  differentiable  in  the  feasible  region. 

For  this  study  f(x)  =  0  at  all  x  and  therefore  none  of  the  nonlinear 
subroutines  were  used  for  problem  solution. 

For  linear  programs,  MINOS  uses  the  revised  simplex  algorithm 
with  all  data  and  instructions  residing  in  core  storage.  The  basis 
inverse  is  maintained  as  an  LU  factorization  using  a  Bartels-Golub  update. 
The  reinversion  routine  uses  the  Hellerman-Rarick  [11]  pivot  agenda 

algorithm.  _ _  _ 

XMP  is  a  library  of  FORTRAN  subroutines  which  can  be  used  to  solve 
linear  programs.  The  basis  inverse  is  maintained  in  LU  factored  form. 

The  pricing  routine  uses  a  candidate  list  of  size  6  with  two  hundred 
columns  being  scanned  each  time  the  list  is  refreshed.  The  basis  is 
reinverted  every  50  iterations. 
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LISS  stands  for  "Linear  In-Core  Simplex  System"  and  is  an  in-core 
LP  solver  with  the  basis  inverse  maintained  in  product  form.  The 
reinversion  routine  is  a  modification  of  the  work  of  Hellerman  and 
Rarick  [11].  The  basis  inverse  is  refactored  every  50  iterations.  A 
partial  pricing  scheme  is  used  with  20  blocks. 

MCNF  stands  for  "Multicommodity  Network  Flow".  .:CNF  uses  the  primal 
partitioning  algorithm  also.  The  basis  inverse  is  maintained  as  a  set  of 
rooted  spanning  trees  (one  for  each  commodity)  and  a  working  basis  inverse  in 
product  form.  This  working  basis  inverse  has  dimension  equal  to  the  number 
of  binding  GUB  constraints.  A  partial  pricing  scheme  is  used.  Our  computa¬ 
tional  experience  is  given  in  Table  1. 

The  row  entitled  GUB  Constraints,  gives  the  number  of  LP  rows  which 
correspond  to  "GUB  Constraints".  The  row,  entitled  "Binding  GUB 
Constraints",  gives  the  number  of  GUB  constraints  met  as  equalities  at 
optimality  using  MCNF .  All  runs  were  made  on  the  CDC  6600  at  Southern 
Methodist  University  using  the  FTN  compiler  with  the  optimization  feature 
enabled. 

Based  on  these  results,  we  conclude  that  for  lightly  constrained 
multicommodity  network  flow  problems 

(i)  XMP  and  MINOS  run  at  approximately  the  same  speed, 

(ii)  NETSIDE1,  NETSIDE2  and  NETSIDE3  run  at  approximately  the  same 
speed,  and 

(iii)  the  three  NETSID  codes  are  approximately  twice  as  fast  as  XMP 


and  MINOS. 
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ABSTRACT 


This  paper  describes  a  computational  experiment  comparing  a 
pure  network  code,  IBM's  MPSX/370,  and  our  implementation  of  a 
heuristic  version  of  the  projective  transformation  algorithm  first 
suggested  by  N.  Karmarkar.  On  five  randomly  generated  dense  assign¬ 
ment  problems,  we  found  that  the  pure  network  code  was  18  times 
faster  than  MPSX  which  was  14  times  faster  than  our  projective 
tranforraation  code. 
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I .  INTRODUCTION 


This  paper  describes  a  pure  network  implementation  of  the  new 
projective  transformation  algorithm  for  linear  programs  developed  by 
N.  Karmarkar  [2],  The  projection  of  the  gradient  in  the  transformed 
space  is  accomplished  by  solving  a  least  squares  problem  using  the 
LSQR  routine  of  Paige  and  Saunders  [4].  On  dense  assignment  problems, 
we  found  that  a  pure  network  code  NETFLO  [3]  is  approximately  250 
times  faster  than  the  new  code  and  MPSX/370  is  14  times  faster  than 
the  new  code.  Other  computational  experience  with  an  early  version  of 
the  algorithm  may  be  found  in  Tomlin  [5], 
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II.  THE  ALGORITHM 


Let  the  linear  program  be  given  by 


min  cx 

(1) 

s.t.  Ax  =  b 

(2) 

£  <  x  <  u  . 

(3) 

The  algorithm  which  we  implemented  may  be  stated  as  follows: 

PROJECTIVE  TRANSFORMATION  ALG 

0.  Initiali2ation 

Let  z*  denote  the  optimal  objective  value  of  (l)-(3)  and  let  x 
be  a  starting  point  such  that  Ax  =  b  and  £  <  x  <  u.  Select  the  step 
size  g  with  0  <  g  <  1. 

1.  Form  Transformation  Matrix 

dj  MIN(u^  -  x ^  ,  x^  -  £  j ) ,  D  =  diagCd^,  d^)  . 

2.  Transform  Constraints 
B  =  AD 

3.  Project  Gradient 

c  «  Dc  -  B"(BB")-1BDc 

4.  Transform  To  Original  Space 
h  =  Dc 

5.  Determine  Max  Step  Size 


0^  -  MIN 

hj  >0 


x .  -  £ . 
J  J 
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a  =  MIN  (o^,  a2 ) 

Move  To  New  Point 
x  =  x  -  a  8  h 

Check  For  Termination 

If  .9cx  <  z* ,  stop  with  x  a  near  optimum 
otherwise,  go  to  1. 


III.  THE  CODE 


We  developed  a  FORTRAN  code,  called  PTANET,  for  the  projective 
transformation  algorithm  which  was  specialized  for  pure  network 
problems.  That  is,  A  is  assumed  to  be  a  node-arc  incidence  matrix, 
less  one  row.  Hence,  A  has  full  row  rank.  Step  3  was  performed  using 
the  subroutine  LSQR  developed  by  Paige  and  Saunders  [4],  That  is,  we 
solve  the  following  least  squares  problem, 

min  |]  B  x  -  Dc  ||  2 

to  obtain  c.  We  used  the  1978  version  of  LSQR  since  that  version 
returns  the  residual.  The  only  two  calculations  involving  B,  p  =  Bv 
and  p  =  u”*B,  were  performed  by  special  routines  which  exploited  the 
network  structure  of  B. 

LSQR  uses  an  input  parameter,  EPS,  for  termination  of  the  least 
squares  solution.  An  EPS  of  l.E-8  was  found  to  be  too  large.  That 
is,  an  x  was  generated  in  which  at  least  one  component  of  |Ax  -  bj 
was  greater  than  l.E-6.  Similar  problems  were  encountered  when  we 
ran  LSQR  in  single  precision.  Hence,  all  arrays  are  double  precision 
and  we  used  the  following  tolerances  and  input  limits  for  LSQR: 

EPS  =  1.0E-12 

ATOL  *  EPS*1000. 

BTOL  -  EPS*1000. 

CONLIM  =  1 . / (10. *DSQRT(EPS) ) 

ITNLIM  *  1000. 
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The  core  storage  comparison  between  PTANET  and  the  pure  network 


code  NETFLO  [ 3 ]  is  as  follows: 


Code 

NETFLO 


Arc  Length  Arrays 


3 


PTANET 


10 


The  1982  version  of  LSQR  requires  fewer  arrays.  The  above  implementation 
of  PTANET  does  not  have  the  minimum  number  of  arrays  that  can  be  achieved. 
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IV.  THE  EXPERIMENT 


Due  to  the  fact  that  our  algorithm  requires  a  starting  point  x 
such  that  Ax  =  b  and  £  <  x  <  u,  we  restricted  our  test  problems  to 
dense  assignment  problems.  That  is,  an  assignment  problem  with  N  =  10 
has  20  nodes  and  100  arcs.  The  costs  were  randomly  generated  integers 
on  the  interval  (1,  100).  For  a  problem  of  size  N,  the  starting 
solution  was  x^  =  1/N  for  all  j. 

We  solved  all  problems  using  NETFL0  first.  The  optimal  objective 
value  was  then  fed  to  PTANET  to  use  for  termination.  We  ran  PTANET 
with  B  =  0.9.  We  also  ran  PTANET  which  rounded  to  the  nearest  feasible 
solution  whenever  N  arcs  had  flow  of  at  least  0.5. 

The  same  problems  were  also  run  on  MPSX/370  using  the  default 
parameter  settings.  Since  two  runs  on  the  same  problem  may  take  a 
different  number  of  iterations,  we  ran  each  problem  three  times  and 
reported  the  average  of  these  runs. 

Our  results  are  given  in  Table  1.  All  runs  were  made  on  the  IBM 
3081-D24  at  Southern  Methodist  University.  NETFL0  and  PTANET  are 
written  in  FORTRAN  and  were  run  using  FORTVS  with  OPT  =  3.  NETFLO 
solved  all  5  problems  in  less  than  1  second,  MPSX  took  18  seconds, 
while  PTANET  required  255  seconds.  The  final  three  iterations  for  each 
run  with  PTANET  required  approximately  seventy  percent  of  the  total 
computational  time.  This  is  due  to  the  ill-conditioning  of  B  =  AD.  As 
the  flows  approach  their  bounds,  the  components  of  D  become  quite  small. 
At  optimality,  every  flow  is  either  at  its  upper  or  lower  bound. 
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V.  THE  NULL  SPACE  MATRIX 


It  appears  that  the  main  computational  problem  with  the 
Projective  Transformation  ALG,  as  stated  in  Section  II,  is  that 
if  there  is  error  in  the  calculation  of  d  in  Step  4,  then  x 
becomes  infeasible.  That  is,  Ax  #  b.  Therefore,  many  iterations 
are  required  by  LSQR  to  obtain  a  sufficiently  accurate  d  so  that 
feasibility  is  maintained.  In  order  to  overcome  this  numerical 
problem,  we  modified  the  algorithm  to  use  the  null  space  matrix  to 
accomplish  the  projection. 

Recall  that  the  direction  is  obtained  by  the  following  steps: 

2.  B  =  AD. 

3.  £  -  (I-B"(BB>)'1B)Dc. 

4.  d  =  Dc . 

Suppose  A  is  m  x  n  and  let  Q-*  be  the  (n  x  n  -  m)  null  space  matrix 
corresponding  to  A.  The  null  space  matrix  corresponding  to  B  =  AD  is 
D  .  By  the  property  of  the  null  space  matrix, 

i-b'(bb')~1b  =  d~1q>(qd"1d_1q')~1qd_1. 

Hence,  2,  3,  and  4  can  be  replaced  by  the  following: 

2.  Construct  Q". 

3.  c  «  D_1Q(QD_1d“V)_1QD-1Dc. 

4.  d  =  Q'(QD~1D~1Q'')“1QD'1Dc. 

By  applying  LSQR  to 

min  jj  D  ^Q'x  -  Dc  ||2 , 


we  obtain 


x*  *  (QD-1D~1Q,*)“1QD''1Dc. 

Then  d  is  simply  Q"x*. 

For  pure  network  problems,  Q  can  be  generated  from  A  and  any 
basis  (rooted  spanning  tree).  The  columns  of  Q'  may  be  constructed 
by  tracing  cycles  in  the  basis  tree  after  a  nonbasic  arc  is  appended 
to  the  tree.  We  modified  PTANET  to  use  Q  and  developed  special 
routines  to  calculate  D  ^Q"*x  and  y“*QD  ^  required  by  the  1982  version 
of  LSQR. 

Computationally,  we  found  that  an  inaccurate  x*  from  LSQR  still 
produced  a  feasible  direction  d.  However,  these  directions  would 
not  necessarily  guarantee  that  optimality  to  the  original  problem 
could  be  obtained.  Our  experience  with  a  20  x  20  assignment  problem 
is  presented  in  Table  2.  With  a  tolerance  of  l.E-10  and  smaller, 
convergence  was  achieved.  However,  with  a  tolerance  of  l.E-8,  the 
objective  function  stalled  at  an  objective  value  13%  above  the  true 
optimum.  That  is,  we  were  trading  a  worse  direction  in  exchange  for 
a  guaranteed  feasible  direction.  Since  the  dimension  of  Q"  was  very 
large  and  D  ^Q"  becomes  ill-conditioned  after  a  few  major  iterations, 
this  trade-off  does  not  pay  off.  NETFLO  solved  this  problem  in  0.02 
seconds. 
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VI.  CONCLUSIONS 


Based  on  our  computational  tests,  we  were  unable  to  confirm 
the  original  claims  which  were  reported  concerning  this  algorithm. 

The  biggest  difficulty  appears  to  be  that  both  AD  or  D  become 
very  ill-conditioned  as  components  of  the  solution  vector  approach 
either  their  upper  or  lower  bounds.  Two  approaches  have  been 
suggested  to  help  alleviate  this  problem.  When  a  variable  gets 
close  to  either  its  upper  or  lower  bound,  fix  it  to  the  appropriate 
bound  and  drop  it  from  the  problem.  Scaling  should  also  assist  in 
this  difficulty.  We  also  found  that  skipping  steps  2,  3,  and  A  and 
using  the  same  direction  two  successive  major  iterations  can  reduce 
the  computational  time  by  up  to  25%.  However,  to  make  our  present 
code  competitive  with  MPSX/370  these  ideas  would  have  to  perform 
spectacularly. 

The  problem  of  a  feasible  interior  starting  point  can  be  solved 
by  a  two-phase  approach.  The  problem  of  a  satisfactory  stopping  rule 
can  be  solved  by  iterating  between  the  primal  and  the  dual.  When  the 
two  bounds  are  within  a  given  tolerance,  then  the  algorithm  terminates. 
These  procedures  could  increase  the  computational  times  by  a  factor 


of  four. 
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Table  1.  Computational  Experience  With  Dense  Assignment  Problems 

(Arcs  =  NxN,  Nodes  =  2N) 
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The  data  in  these  rows  is  identical;  the  heuristic  was  unable  to  produce  an  early  termination 
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ABSTRACT 


This  paper  gives  a  mathematical  programming  model  for  the  problem 
of  assigning  frequencies  to  nodes  in  a  communic itions  network.  The 
objective  is  to  select  a  frequency  assignment  which  minimizes  both  co¬ 
channel  and  adjacent  channel  interference.  In  addition,  a  design  engineer 
has  the  option  to  designate  key  links  in  which  the  avoidance  of  jamming 
due  to  self-interference  is  given  a  higher  priority.  The  model  has  a 
nonconvex  quadratic  objective  function,  generalize!  upper  bounding  con¬ 
straints,  and  binary  decision  variables.  We  developed  a  special  heuristic 
algorithm  and  software  for  this  model  and  tested  it  on  five  test  problems 
which  were  modifications  of  a  real-world  problem.  Even  though  most  of  the 
test  problems  had  over  600  binary  variables,  we  were  able  to  obtain  a 
near  optimum  in  less  than  12  seconds  of  CPU  time  on  a  CDC  Cyber-875. 
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I.  INTRODUCTION 


One  of  the  most  critical  design  problems  in  a  radio  communication 
network  is  the  assignment  of  transmit  frequencies  to  stations  (nodes)  so 
that  designated  key  communication  links  will  not  be  jammed  due  to  self¬ 
interference.  In  this  investigation,  we  describe  a  novel  new  optimization 
model  and  a  solution  technique  which  can  be  used  to  assist  design  engineers 
in  this  process. 

1.1.  Problem  Description 

A  radio  communications  network  consists  of  radio  stations,  each 
equipped  with  one  or  more  transmitters  and  receivers.  When  a  given 
station  has  the  ability  to  receive  Information  intelligibly  from  a  trans¬ 
mitting  station,  a  link  is  said  to  exist  from  the  transmitting  station  to 
the  receiving  station.  The  interconnection  of  these  stations  and  links 
may  be  viewed  graphically  as  a  set  of  nodes,  representing  the  radio  stations, 
joined  together  by  directed  arcs,  representing  the  links. 

We  assume  in  our  model  that  one  transmitter  and  several  receivers  are 
locatd  at  each  radio  station  (node).  The  transmitter  Is  tuned  to  a 
specified  center  frequency,  and  the  receivers  are  tuned  to  the  transmit 
frequencies  of  the  neighboring  stations  to  which  the  station  is  to  be 
linked.  A  channel  is  associated  with  each  center  frequency  in  a  way 
similar  to  the  way  channels  and  frequencies  are  associated  in  a  television 
set.  When  a  TV  is  tuned  to  channel  4,  for  example,  it  is  really  being 
tuned  to  receive  video  signals  being  broadcast  at  67.25  MHz. 
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For  our  model,  a  given  center  frequency  will  be  associated  with 
each  channel  number.  Using  this  definition,  the  frequency  assignment 
problem  may  be  defined  as  follows: 

"Given  N  transmitting  stations  (nodes),  assign  1  of  F  transmit 
channels  to  each  node  in  such  a  way  as  to  minimize  the  number 
of  designated  key  links  jammed  due  to  co-channel  and  adjacent 
channel  interference," 

We  say  that  a  link  is  jammed  if  either  of  the  following  conditions 
occurs : 

(i)  a  node  receives  two  signals  on  the  same  channel  that  are  less 
than  a  dB  apart  in  signal  strength,  or 

(ii)  a  node  receives  a  signal  on  a  given  channel  while  a  neighboring 
node  transmits  on  an  adjacent  channel.  If  the  neighbor's 
signal  strength  exceeds  the  signal  strength  of  the  current  node 
by  more  than  3  dB,  then  the  incoming  signal  will  be  garbled. 

The  constants  a  and  3  are  functions  of  the  hardware  used  in  the 
network.  Some  of  the  determining  factors  are  the  receiver  selectivity, 
the  type  of  signal  modulation,  and  the  purity  of  the  signal. 

We  now  introduce  the  notation  used  to  describe  the  mathematical 
model.  Let  f  c  {1,  ...,  F}  denote  a  channel  and  n  e  {1,  ...,  N}  denote 
a  node,  e^  will  denote  a  vector  whose  entries  are  0  except  for  the  ic^ 
which  is  1.  Let 

x,  ■  1  if  channel  f  is  assigned  to  node  n 

fn 

and  0  otherwise. 
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-  the  row  vector  [x^,  •••»  x£ »  and 

g(x. ,  ,,,,  x„)  ■  a  weighted  number  of  jammed  links  with 
—1  — F 

assignment  (x^,  . ..,  Xp) • 

Using  the  above  notation,  the  mathematical  model  of  the  frequency 
assignment  problem  is 


min 

g(xx,  •••.  XF> 

(1) 

s.t. 

Ex,  »  1,  all  n 
f  fn 

(2) 

Xfn  £  (0,1),  all  f,n. 

(3) 

For  this  application,  g(*)  is  a  nonconvex  quadratic  function  and 
therefore  (1)  -  (3)  is  a  member  of  the  class  of  binary  nonconvex  cost 
nonlinear  programs. 

1.2  Related  Literature 

A  heuristic  procedure  for  solving  a  similar  problem  using  a  graph 
coloring  algorithm  has  been  evaluated  by  Zoellner  and  Beall  [7].  Closely 

related  models  have  been  investigated  by  Morito,  Salkin,  and  Williams 
[5)  and  by  Mathur,  Salkin,  Nishimura,  and  Morito  [4],  Their  models 
are  general  linear  integer  programs  with  a  single  constraint.  Using 
a  special  branch-and-bound  algorithm,  they  successfully  solved  their 
model  with  up  to  fifty  channels. 

1.3  Accomplishments  of  the  Investigation 

We  developed  a  novel  new  mathematical  model  of  the  frequency 
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assignment  problem  which  takes  the  form  of  a  binary  nonconvex  quadratic 
cost  nonlinear  program.  The  model  incorporates  weighting  constants  that 
allow  a  design  engineer  to  tune  the  model  to  a  particular  application. 

We  present  an  elegant  specialization  of  the  convex  simplex  algorithm 
to  obtain  a  local  optimum  for  this  model.  In  addition,  specialized 
software  has  been  developed  for  this  model  and  tested  on  five  versions 
of  a  real-world  problem.  The  software  works  quite  well  requiring  less 
than  a  minute  of  computer  time  for  all  five  test  problems. 


4 


II.  THE  OBJECTIVE  FUNCTION 


In  this  section,  we  define  the  weighted  interference  function, 
g(x  ,  ...,  x  ).  This  function  is  generated  from  a  set  of  signal 
strength  matrices,  (A^,  ....  A^,) ,  two  weighting  matrices,  and  a  set 
of  critical  values  a,  $,  and  6^,  ...»  <5^.  Let  a^  denote  the 
received  signal  strength  in  dBu/m  of  a  signal  which  originates  at 
node  i  and  is  received  by  node  j,  and  let  A^  denote  the  matrix  whose 
elements  are  a^.  Let  the  weighting  matrices  P  and  W  be  determined 
as  follows: 


ij 


and 


ij 


p^,  if  (i,j)  is  a  designated  key  link 
P2,  otherwise 

w^,  if  (i,j)  is  a  designated  key  link 


w2>  otherwise. 


The  constants  p^,  p^,  w^,  and  w^  are  tuning  parameters  which  are  used 
to  provide  weights  in  the  interference  function  for  the  key  links. 
Gamma  is  used  to  denote  the  scalar  function,  defined  by 

1,  if  x  >  0 

1 

Y(x)  = 

) 

0,  otherwise. 
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Using  y(«)»  we  define  the  three  matrices 


’  k  j  -  4|,wik +  1  * i 

aL. 5  \ 


0,  otherwise; 


E  Y(a^  '  4 ' 6)  Wik»  1  *  J 


k  +  i,j 

4  >  \ 


0,  otherwise; 


'  „  1 1 ,  -  4  - B)  "ik> 1  * J 

k  t  i»j 

!a!.  >  <5, 

ik  k 


0,  otherwise. 

Using  these  matrices,  the  interference  function  is  given  by 

g(xr  ...»  Xp) 


Z  x'  Qf  xf 
f=l  " 


f-F-1 

Z  x'  Rf  xf+1 
f=l 


f9  -f  Sf  -f-1  ' 

f*2 


co -channel 
interference 


adjacent  channel 
interference  from 
channel  above 


adjacent  channel 
interference  from 
channel  below 


In  addition  it  is  often  desirable  to  use  all  of  the  channels. 
Therefore,  we  appended  the  function 


i=N-l  j»N 


w  Z  I  x  I  x 
J  f  i*l  1  j=i+l  13 


to  g(«)  so  that  in  the  absence  of  self-interference,  the  channels 
would  be  equally  distributed  among  the  nodes.  The  scalar  w^  is  also 
a  tuning  parameter. 

Using  the  above  formulae,  we  now  give  an  example  which  presents 
the  matrices  required  to  define  g(»).  Let  a  **  2,  B  *  3,  **  0,  *  0 

for  all  n,  and  *  w^.  *  1  for  all  i,j. 


0  12  5 


10  3  3 


2  3  0  2 


4  3  2  OJ  ,  and 


0  2  3  5 


2  0  5  5 


3  5  0  1 


5  5  1  Ol  ,  then 


0  10  1 


10  2  1 


0  2  0  1 


1110 


0  10  0 


10  10 


0  10  1 


0  0  10 


0  0  11' 


0  0  0  1 


0  0  0  0 


0  0  0  0  J  ,  and 


0  0  0  0' 


0  0  0  0 


10  0  0 


0  0  0  0 


■*  ■  jVv 


Ill .  THE  ALGORITHM 


Let  x"*  *  [x',  xp.  Then  the  frequency  assignment  problem 

takes  the  general  form: 


min 

g(x)  =  x"  C  x 

(A) 

(0 

• 

ft 

• 

Z  xr  ml,  all  n 
f  fn 

(5) 

Xfn  £  {0,1},  all  f,n 

(6) 

where  the  diagonal  elements  of  C  are  0  and  all  other  elements  are  positive. 

The  continuous  relaxation  of  (A)  -  (6)  is  obtained  by  replacing  (6)  with 

°  <  x^  <  1,  all  f  *n.  (7) 

The  model  (A),  (5),  (7)  is  a  nonconvex  quadratic  program  and  a  local 

optimum  can  be  efficiently  obtained  by  application  of  the  convex  simplex 

algorithm  as  described  in  Zangwill  [6j.  Suppose  we  begin  with  a  feasible 

integer  solution  x"  *»  [x^,  xp.  We  assume  that  all  nonbasic  variables 

have  a  value  of  zero.  Let  ...,  denote  the  subscript  such  that 

x„,«...=x  •  1.  Then  a  nonbasic  variable  x,  with  a  value  of  zero, 

JL1  l  n  fn 

1  n  _ 

prices  favorably  if  [Vg(x)]'*  (e^  -  ep  <  0  where  i  *  (f  -  1)  N  +  n  and 
j  *  (f  -  1)  N  +  The  line  search  for  this  problem  requires  that  we 

solve  the  problem 

min  g(x  +  (e.  -  e ,)A).  (8) 

0  <  A  <  1  " 

But  dg(x  +  -  e^A) 

dA 

•  (e4  -  ep "  Vg  (x  +  ei  -  ep 

■  <Si  -  £j»'  <c  +  n  <i  +  u  -  «j) 


A  -  1 
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-  (Vg(x)  J "  (e±  -  ej)  +  (et  -  e (C  +  C')  (e±  -  £j). 

Since  xfn  priced  favorably,  then  [Vg(x)K  (e^  -  e^)  <  0. 

Also,  (e^  -  e^)  "  (C  +  C'*)  (e^  -  e^) 

-  e'  (C  +  C')  e±  +  e'  (C  +  C')  e  -  e'  (C  +  C')  £j  -  e'  (C  +  C')  e^ 
But,  the  diagonal  elements  of  (C  +  C')  are  0  and  all  other  elements  are 
nonnegative.  Hence,  the  solution  to  (8)  is  4*  ■  1  and  the  exact  change 
to  the  objective  function  will  be  Vg  (x')  (e^  -  e^)  -  e^"  (C  +  C")  e^ 

-  ej"  (C  +  C *)  e^  ,  a  strict  decrease.  Therefore,  in  the  new  solution 

x,  is  set  to  1  and  x.  is  set  to  0.  Since  this  holds  for  every 
fn  i-  n  J 

n 

iteration  of  the  convex  simplex  algorithm.  Integrality  is  maintained 
and  a  local  optimum  for  (4)  -  (6)  can  be  obtained  by  finding  a  local 
optimum  for  (4),  (5),  (7). 

Let  x  be  any  initial  assignment  for  the  frequency  assignment 

problem.  Using  this  Initial  assignment,  the  algorithm  may  be  stated 

as  follows: 

For  f  -  1,  ....  F. 

For  n  »  1,  ...,  N. 

L  :  •  k  where  x.  ■  1. 
n  kn 

i  s  -  (f  -  1)  N  +  n. 

j  :  -  (f  -  1)  N  +  £n. 

P  *  «  [Vg(x)r  <£l  -  £j). 

If  p  <  0 

then 


Repeat  as  long  as  p  <  0  for  some  f  and  some  n. 
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IV.  COMPUTATIONAL  EXPERIENCE 


We  implemented  the  frequency  assignment  algorithm  in  a  FORTRAN 
code.  All  data,  including  the  matrices  Q^,  R^,  and  S^,  are  stored  in 
high  speed  core.  Special  subroutines  were  written  to  evaluate  both 
g(«)  and  Vg(*)  at  a  point.  The  code  begins  with  F  different  starting 
solutions  and  stops  when  a  local  optimum  is  found.  The  initial  assign¬ 
ment  for  run  r  £  {l,  ...,  F}  is  to  assign  frequency  {[(n  +  r  -  2)  modulo  F] 

+  l}  to  node  n.  The  best  solution  obtained  from  all  F  runs  is  the 

output. 

Five  test  problems  were  generated  from  the  real-world  43  node  network 
illustrated  in  Figure  1.  The  lines  connecting  nodes  are  the  designated 
key  links.  The  problems  all  have  the  same  topology  but  differ  in  the 
selection  of  the  critical  values  and  the  weighting  constants.  A  random 
assignment  was  generated  and  the  matrices  were  modified  so  that  this 
assignment  produced  a  cost  of  zero.  Hence,  the  optimal  objective  value 
for  each  problem  is  zero. 

Our  computational  experience  is  reported  in  Table  1.  All  runs 
were  made  on  a  CDC  Cyber-875  using  the  FTN5  compiler  with  OPT  **  2. 

The  "Initial  Obj  Value"  row  is  the  average  objective  value  for  the  F 

initial  solutions.  Note  that  all  five  problems  were  run  in  less  than 
1  minute  of  CPU  time  and  the  "Final  Obj  Values"  were  quite  close  to  the 
optimum  as  compared  to  the  initial  assignments. 


Figure  1  Table  1 
About  Here 


V.  CONCLUSIONS 


Our  optimization  model  and  computer  software  provide  a  practical 
approach  to  assist  communication  network  designers  in  obtaining  near 
optimal  solutions  for  the  frequency  assignment  problem.  The  fact  that 
the  diagonal  elements  of  C  in  the  quadratic  objective  function  x^  C  x 
are  zero,  allows  a  very  efficient  implementation  of  the  convex  simplex 
method  which  maintains  integrality.  Hence,  if  we  begin  with  an  integer 
assignment,  the  convex  simplex  algorithm  follows  a  sequence  of  integer 
points  until  a  local  minimum  is  obtained.  This  procedure  is  so  fast 
that  very  large  problems  can  be  easily  handled  by  this  approach. 
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Table  1.  Computational  Results  With  43  Node  Model 


Row 

Description 

Problem 

1 

2 

3 

4 

5 

ctdB 

10 

10 

12 

10 

10 

BdB 

25 

25 

25 

25 

30 

F  (channels) 

10 

12 

14 

14 

14 

Binary  Variables 

430 

516 

602 

602 

602 

Iterations 

525 

341 

497 

540 

526 

Solution  Time  (secs) 

5 

7 

10 

11 

11 

Initial  Obj  Value 

3153 

1561 

1875 

1671 

1670 

Final  Obj  Value 

164 

103 

79 

5 

4 

Jammed  Key  Links 

8 

4 

_ 

3 

0 

0 

Figure  1.  Forty-three  Node  Communication  Network  With  Designated  Key  Links 


